# How to plot the contour of f[x,y]==0 if always f[x,y]>=0

I've been stuck at this problem for weeks and I've asked several related questions here:

ContourPlot shows only part of the contours [duplicate]

How to plot the contour of the radius part of a complex function

a lot of help has been kindly given but I think my problem is still not solved completely.

In the first question, Jens pointed out that when a function only touch but not cross 0, ContourPlot would have difficulty finding the 0 contour. He then suggested to plot f[x,y]==0.001 or some small number, instead of f[x,y]==0.. This seems works very well for some functions like Abs[Cos[x]+Cos[y]]. However it doesn't work on my function, and he also gave an explanation that the reason it doesn't work is because the maximals of my function change a lot in the plot region. Along this direction, I rescaled my function to make the maximals at the same order, but the f[x,y]==0.001 trick still doesn't work. I've been experimenting with these for a while and found that the reason seem more than just the steep change in the maximals. And I also found some behavior of ContourPlot that are very confusing for me.

Consider this function, the maximals of it changes quite of lot

f[x_, y_] := (Cos[x] + Cos[y]) Exp[I x y] Exp[x/2]

Plot3D[Abs[f[x, y]] == 0.001, {x, 0, 4 Pi}, {y, 0, 4 Pi},
PlotPoints -> 40, PlotRange -> All]


Now if we plot the 0 contour, using different PlotPoints setting, we would get different results. In fact more plot points makes the contour worse. Why?

myPlot[sqs__] := Module[{},
Row@{ContourPlot[Sequence[sqs], ImageSize -> 400],
ListPlot[
Reap[ContourPlot[Sequence[sqs],
EvaluationMonitor :> Sow[{x, y}]]][[-1, 1]],
PlotStyle -> PointSize[Tiny], ImageSize -> 400, Axes -> False,
Frame -> True, AspectRatio -> 1]}
]

myPlot[Abs[f[x, y]] == 0.001, {x, 0, 4 Pi}, {y, 0, 4 Pi}]
myPlot[Abs[f[x, y]] == 0.001, {x, 0, 4 Pi}, {y, 0, 4 Pi}, PlotPoints -> 50]
myPlot[Abs[f[x, y]] == 0.001, {x, 0, 4 Pi}, {y, 0, 4 Pi}, PlotPoints -> 100]


If we increase the plot range, the ContourPlot just fails to give the result. Why? The MaxRecursion seems can help in this case, but the contour is still broken.

myPlot[Abs[f[x, y]] == 0.001, {x, 0, 16 Pi}, {y, 0, 16 Pi}]
myPlot[Abs[f[x, y]] == 0.001, {x, 0, 16 Pi}, {y, 0, 16 Pi}, MaxRecursion -> 3]


Now compare with the second function

f[a_, a0_, k_, K0_] :=
a^2 Sech[(a a0)/2]^2 (-2 I (1 + E^(2 I a0 k)) k +
a (-1 + E^(2 I a0 k)) Tanh[(a a0)/2]) +
2 k (I E^(
I a0 k) ((a - k) (a + k) Cos[a0 k] + (a^2 + k^2) Cos[a0 K0]) +
a (-1 + E^(2 I a0 k)) k Tanh[(a a0)/2])

With[{a0 = 10., a = 1.4},
Plot3D[Abs[f[a, a0, k, K0]] == 0.001, {K0, -2 π/a0,
2 π/a0}, {k, 0, 2}, ImageSize -> 400, PlotPoints -> 40]]


We can see that this function behaves well in terms of the changes in the maximals. All the maximals seems in the same order. And ContourPlot seems fails to produce the results, even plotting a very small region with MaxRecursion->3. Why?

With[{a0 = 10., a = 1.4},
myPlot[Abs[f[a, a0, y, x]] == 0.001, {x, -π/a0, π/a0}, {y, 0,2}]]
With[{a0 = 10., a = 1.4},
myPlot[Abs[f[a, a0, y, x]] == 0.001, {x, -π/a0, π/a0}, {y, 0,2}, MaxRecursion -> 3]]


Summary of the Questions

1. In the first function, why increase PlotPoints makes things worse, and why increases plot range resulted in almost no contour plotted(when plot range is 16pi)?
2. Why the second function, which seems behave better than the first, shows very poor contour results and how to improve it?
-
Related but not the same question: ContourPlot shows only part of the contours –  ybeltukov Sep 20 '13 at 21:50
–  Alexey Popkov Sep 26 '13 at 7:57

My other answer addressed the question in the title. Here I will try to address the two queries at the end of the OP question. @ybeltukov was the first to point out the key issue, which is that ContourPlot works only when it finds sample points where the values of the function are greater than the contour level and points where they are less than it.

ContourPlot, like other plot functions, starts with a 2D rectangular grid of sample points determined by PlotPoints, which is then recursively subdivided up to MaxRecursion times. Where or when a recursive subdivision occurs depends on the function and the contour levels. Some discussion of how the subdivision works can be found in the answers to Specific initial sample points for 3D plots.

I hope to show that the answer to the first question about why increasing PlotPoints results in a worse graph is mainly due to the alignment of the sample points with a certain region, the extremely small region where the function is less than 0.001. The region is so small that misalignment is possible even when PlotPoints is set well over 100.

### Auxiliary functions

There are a couple of helper functions (cpGrid, cpShow) below that create ContourPlots of a function that show how the plot domain is subdivided together with the sample points at which the values of the function are less than the contour level. (See "Code dump" below at the end.)

The OP's functions, which I will call f1 and f2:

ClearAll[f1, f2];
SetAttributes[f1, Listable];
SetAttributes[f2, Listable];
f1[x_, y_] := (Cos[x] + Cos[y]) Exp[I x y] Exp[x/2];
f2[a_, a0_, k_, K0_] :=
a^2 Sech[(a a0)/2]^2 (-2 I (1 + E^(2 I a0 k)) k +
a (-1 + E^(2 I a0 k)) Tanh[(a a0)/2]) +
2 k (I E^(I a0 k) ((a - k) (a + k) Cos[a0 k] + (a^2 + k^2) Cos[
a0 K0]) + a (-1 + E^(2 I a0 k)) k Tanh[(a a0)/2]);


## ContourPlot sampling

The table of plots below show PlotPoints settings of 3, 4, 5 points (rows) and MaxRecursion settings of 0, 1, 2 (columns). The initial sample points lie on a rectangular grid, but the actual subdivision is made of triangles with a SW-NE bias. This is shown below by the gray Mesh lines. The red points are sample points where the value of the function Abs@f1 lies below the level 0.001. The other intersections are sample points where the value of the function is above the level. Where a red point is next to a gray intersection, one may observe an increase in the number of subdivisions.

{cp, samp} = cpGrid[Abs@f1@## &, {x, 0, 4 Pi}, {y, 0, 4 Pi}, 0.001, {3, 4, 5}];
cpShow[cp, samp, Abs@f1@## &, 0.001] // GraphicsGrid


There are several things to observe. The zeroes of f1 lies on lines where $x \pm y$ is a multiple of $\pi$. For the domain $[0,4\pi]\times[0,4\pi]$, the sample points will happen to lie on these lines when PlotPoints is of the form $4n+1$. So we see red points for 5 points with no subdivision, but not initially in the other plots. Below are some more examples showing the contour plots with the initial sample points (MaxRecursion -> 0) that exhibit the $4n+1$ pattern. Compare 100 vs. 101 and 25 vs. 101.

GraphicsGrid@Table[
With[{tolerance = 0.001, pp = p + 4 j},
Show[ContourPlot[
Abs[f1[x, y]] == tolerance , {x, 0, 4 Pi}, {y, 0, 4 Pi},
PlotPoints -> pp, MaxRecursion -> 0, PlotRange -> All,
PlotLabel -> pp], PlotRange -> {{0, 4 Pi}, {0, 4 Pi}}]
],
{j, {1, 2, 20}}, {p, 20, 23}
]


Now a single subdivision of every triangle would transform a grid of $n$ by $n$ points into one with $2n-1$ plot points in each direction, which, if repeated, would eventually bring all grids into the form of $4n+1$ points on a side; however depending on the function and contour level, some triangles are subdivided twice and some only once. For example, there is a pair at approximately {5, 8} in the plot for 3 points, 1 subdivision, whose shared hypotenuse was not divided even though the midpoint is a zero of the function. I do not know how ContourPlot decides whether to subdivide, but the lack of a division here makes a gap in the contour. (One can check that the contour does not connect these points. A mesh line connecting two red points seems to be interpreted as being connected by a region in which the value of the function is less than 0.001.) In the next step where MaxRecursion increases from 1 to 2, it is subdivided and the gap is closed. One observes a similar gap in the plot for 4 points, 2 subdivisions around the coordinates {8.5, 5}. That gap disappears with another increase to MaxRecursion -> 3. In fact, I think that whatever the number of initial PlotPoints, the contour will be drawn completely with MaxRecursion -> 3 for the symmetric domain $[0,4\pi]\times[0,4\pi]$, but it depends precisely because the special symmetry leads to alignment of the sample points with the zeroes of f1. For an odd number of initial PlotPoints, one only needs MaxRecursion -> 2 to get the whole contour, and for a number of the form $4n+1$, no recursion is necessary.

The dependence on the symmetry can be seen by perturbing the domain slightly:

ContourPlot[
Abs[f1[x, y]] == 0.001, {x, 0, 4 \[Pi] + #}, {y, 0, 4 \[Pi] + #},
MaxRecursion -> 3] & /@ {0., 0.001} // GraphicsRow


When the domain is extended to 16 Pi, the number of plot points needs to be of the form $16n+1$. The following is perhaps the simplest way to get the complete contour (MaxRecursion -> 1 is needed to overcome the SW-NE bias and connect the NW-SE contours):

ContourPlot[Abs@f1[x, y] == 0.001, {x, 0, 16 Pi}, {y, 0, 16 Pi},
PlotPoints -> 17, MaxRecursion -> 1]


Thus the alignment of the grid with the region in $xy$ plane where Abs[f1[x, y]] < 0.001 is the key to the somewhat odd dependence on the number of PlotPoints.

## The thinness of the region in which sample points must land

A single contour level divides the plane into two regions, one over which the value of the function is greater than the level and one over which the value of the function is less, in addition to the level set itself. Normally, if the function is not locally constant, the level set will be a curve. In the OP's examples, we have another sort of edge case in which desired contour curve is the minimum of the function.

When one of the two regions is extremely small, it is unlikely that many sample points will fall into it except "by accident" so to speak. This is what is happening with both of the OP's functions.

To give an illustration we can see, let's use a sequence of values for the level, say, 1., 0.1, 0.01. Despite the images, the true regions are connected, but the width of the region was too narrow to contain enough (or any) points to define a boundary.

GraphicsRow @ RegionPlot[
Abs[f1[x, y]] <= #, {x, 0, 4 Pi}, {y, 0, 4 Pi},
PlotPoints -> 40, MaxRecursion -> 3] & /@ {1, 0.1, 0.01}]


With[{a0 = 10, a = 1.4},
GraphicsRow @ RegionPlot[Abs[f2[a, a0, y, x]] <= #,
{x, -2 Pi/a0, 2 Pi/a0}, {y, 0, 2},
PlotPoints -> 40, MaxRecursion -> 3] & /@ {1, 0.1, 0.01}]]


For the second function, I think it proves to be not as well-behaved as at first glance. It gets quite steep in half of the domain, which accounts for the narrowing of the region as the y coordinate increases. (Note: My plot turns out to have a somewhat greater PlotRange than that of the OP.)

With[{a0 = 10, a = 1.4},
Plot3D[Abs[f2[a, a0, y, x]], {x, -2 \[Pi]/a0, 2 \[Pi]/a0}, {y, 0, 2}, PlotPoints -> 100]]


)

## Summary

I think the behavior of ContourPlot can be accounted in both cases by the extremely small area in the plot domain for which $f < 0.001$. In the case of the OP's first function, there is accidental and sporadic alignment of the zero set of $f1$ with the sample points, which accounts for why a plot with much fewer PlotPoints can look better than one with more.

## Code dump

For investigating contour plots.

cpGrid[f_, xdom_, ydom_, level_, pp0_List, mr0_: 2] :=
Module[{cp, samp},
cp = samp = Table[{}, {pp, Length@pp0}, {mr, 0, mr0}];
Do[{cp[[pp, 1 + mr]] = First@#,
samp[[pp, 1 + mr]] = Hold @@ Last@#} &@
Reap @ ContourPlot[
f[x, y],
{x, xdom[[-2]], xdom[[-1]]}, {y, ydom[[-2]], ydom[[-1]]},
Contours -> {level}, PlotPoints -> pp0[[pp]],
MaxRecursion -> mr, Mesh -> All, EvaluationMonitor :> Sow[{x, y}],
PlotLabel -> Row[{pp0[[pp]], " points, ", mr, " subdivisions"}]],
{pp, Length@pp0}, {mr, 0, mr0}];
{cp, samp}];

ClearAll[cpShow];
SetAttributes[cpShow, Listable];
cpShow[cp_, Hold[samp_], f_, level_] :=
Show[cp,
Graphics[{Red, PointSize[0.02], Point[Pick[samp, UnitStep[f @@@ samp - level], 0]]}]
]

-

Update

I disagree with other answers that blame the problems exclusively on thin $0<f<ε$ or the no touching theory. A counterexample was already given in the last graphs of this answer where the region of interest is about 30% of the plot region. Unless we have access to ContourPlot code, we are really guessing.

In any case, if you want to obtain the following graphs for your first and second examples, follow the recipe I wrote in one of the comments to this answer.

For the first example, I used:

Module[{lf = Log@ComplexExpand@Abs@f1[x, y], xp}, xp = Simplify[D[lf, x]^2 + D[lf, y]^2];
ContourPlot[xp, {x, 0, 4 π}, {y, 0, 4 π}]]


I think that the problem lies in the fact that the regions $0\leq f\leq ε$ have extremely different sizes. Let us analyze the following graph:

myPlot[expression_, contours_: 10] := With[{vect = {D[-expression, y], D[expression, x]}}
, Show[{ContourPlot[expression, {x, 0, 4 Pi}, {y, 0, 4 Pi}
, Contours -> contours, ContourStyle -> None],
StreamPlot[vect, {x, 0, 4 Pi}, {y, 0, 4 Pi},
StreamPoints -> Fine]}]]

myPlot[E^(x/2) Sqrt[(Cos[x] + Cos[y])^2]]


yields

At the right side of the graph, the plot has a lot of contrast and the regions near $f=0$ are rather small. As you go to the left, the graph is flatter and the regions near $f=0$ are rather large. A method that plots contours (including $f=0$) should not solve the equations for the equipotential lines, but instead it should find a value and then "trace its value" within a vertical bin. The vertical size of the vertical bins is determined by the global maximum and minimum; and, that size might not be small enough to discriminate the structure at the left of the graph. Thus, the quality of your graph is not limited by the number of sampling points but by the number of vertical bins. Increasing from 10 vertical bins to 80, one gets a slightly better graph:

myPlot[E^(x/2) Sqrt[(Cos[x] + Cos[y])^2], 80]


Given that your functions come from the absolute value of complex functions, I would enhance the contrast using the facts (1) $f$ is not negative, (2) you are interested in the $f=0$ contour, and (3) $\log(0)=-∞$:

ContourPlot[ Log[E^(x/2) Sqrt[(Cos[x] + Cos[y])^2]], {x, 0, 4 Pi}, {y, 0, 4 Pi}]


Note that the $f=0$ becomes $\log(0)=-∞$ and thus it will be out of any scale. In the previous graph, it shows as absolute white. Since that white is not much different than the other shades of white (which are actually far from the answer), let us change the palette. For your second example, the code

With[{a0 = 10., a = 1.4}, ContourPlot[Log[Abs[f[a, a0, k, K0]]],
{K0, -2 \[Pi]/a0, 2 \[Pi]/a0}, {k, 0, 2},
ImageSize -> 400, PlotPoints -> 80, PlotLegends -> Automatic, ColorFunction ->
"DarkRainbow", ClippingStyle -> Black, Contours -> 80]]


returns

Summary

(1) The quality of a ContourPlot depends not only on the number of plot points but also on the number of vertical bins. The latter is the main problem for the functions you have.

(2) Logarithmic scaling enhances the contrast around $f=0$. If you are interested in the set $f(x)=4$, plot Log[Abs[f[x]-4]].

ContourPlot[f[x,y]==0,blah] does not extract the line f[x,y]==0. It finds an instance of f[x,y]==0 and then it follows it to within a given error: $-ε<f<ε$. The exact algorithm to "follow" it (or "trace" it) is proprietary and we do not know. But we can say some general things. First, notice that in regions where the function is flatter, the algorithm chokes. The effect of the horizontal size of the vertical bins is shown with this code:

Row[{
ContourPlot[If[Abs[x^2 + y^2 - 4] < 1, 1, Abs[x^2 + y^2 - 4]],
{x, -m, m}, {y, -m, m}, ImageSize -> 200],
ContourPlot[If[Abs[x^2 + y^2 - 4] < 1, 1, Abs[x^2 + y^2 - 4]] == 1,
{x, -m, m}, {y, -m, m}, ImageSize -> 200]
}]


For m=2.35, you get a somewhat acceptable graph:

But for m=3, you see that the algorithm begins to wander:

That is, when the horizontal size of the vertical bin is smaller, it is more difficult to find the next point to trace the contour. Then, the algorithm begins to fail.

The above answers your question in the comments: it is not about $ε$ in $0<f<ε$. It is about the size of the neighborhood of $x$ such that $0<f(x)<ε$ (the horizontal size).

-
Thanks for the answer, but 1) If the reason ContourPlot fails is $\epsilon$ changes too much, why Abs[f[x, y]] == 0.001 trick works for the first function but not the second function? We can see from the 3d plots that the first function changes a lot in the $\epsilon$ value while doesn't change much for the second function. 2) Is it possible to extract only the line of f[x,y]==0, since this contour is a dispersion relationship and only makes sense at f[x,y]==0 line. –  xslittlegrass Sep 24 '13 at 14:37
I'm sorry, where did you edit? It shows me your answer was last edited 5 hours ago. –  xslittlegrass Sep 24 '13 at 17:27
By the way, the main reason the solution ContourPlot[Log[Abs[f[x,y]]],blah] works is not because it increases the horizontal size of the vertical bin (which it does) but because the bins accumulate around the singularities. Hence, if you really want to "extract" the lines, you might want to find the places where the gradient of Log[Abs[f[x,y]]]is large. –  Hector Sep 24 '13 at 18:23
I'm confused by your remark about 30% of the plot region the last graphs. In none of the plot region is $f < 1$, where $1$ is the contour level (what I was calling $\epsilon$, which was 0.001 in the OP's question). –  Michael E2 Sep 26 '13 at 2:34

## Introduction

The problem is that the regions where $0 \le f < \epsilon$ can be extremely thin for some functions. ContourPlot works by finding values of $f$ on either side of a contour and refines the contour by recursive subdivision. When ContourPlot doesn't find values on either side of a contour level, it will fail to find the contour curve.

One might hope to detect zeros from the partial derivatives, or as I will use, differences. This is because the slope changes sign at f == 0, if f >= 0. A problem arises if the zero set lies parallel to the step direction, since the slope will be identically zero. In this case one can adjust the step direction. It is always possible to construct a function that will cause trouble, but normally I think most functions can be handled.

## OP's functions

ClearAll[f1, f2];
SetAttributes[f1, Listable];
SetAttributes[f2, Listable];
f1[x_, y_] := (Cos[x] + Cos[y]) Exp[I x y] Exp[x/2];
f2[a_, a0_, k_, K0_] :=
a^2 Sech[(a a0)/2]^2 (-2 I (1 + E^(2 I a0 k)) k +
a (-1 + E^(2 I a0 k)) Tanh[(a a0)/2]) +
2 k (I E^(I a0 k) ((a - k) (a + k) Cos[a0 k] + (a^2 + k^2) Cos[a0 K0]) +
a (-1 + E^(2 I a0 k)) k Tanh[(a a0)/2]);


## Plotting the zero set

We use ContourPlot to get contours where the finite difference f[p + dp] - f[p - dp] changes sign. (Here p = {x, y} and dp is a small step dp = {dx, dy} in some direction. This will gives contours where a path in the direction of dp has a local maximum or minimum. We have to pick out those which are (approximately) along the minimum curve f == 0. We identify which points are approximate zeros as follows. The secant line method is applied to estimate the closest zero of f starting at each given point. The secant line through the given point and a near neighbor is used. If the given point is an approximate zero, the secant line method will produce a point nearby. If the given point is not a near neighbor, it will yield a point far away. (See ptsdf.) The distance from the given point to the zero estimate is used to select the approximate zeros.

The contour lines themselves have to be split into segments that connect approximate zeros and those that do not. (The latter are discarded.) They have to be split because ContourPlot connects lines at random more or less where contour lines intersect.

ClearAll[zeroSet];
Options[zeroSet] =
Join[{"StepSize" -> 0.001, "Tolerance" -> 0.1, "StepVector" -> {0, 1}},
Options[ContourPlot]];
SetAttributes[zeroSet, HoldAll];
zeroSet[f0_, {x_, x1_, x2_}, {y_, y1_, y2_}, opts : OptionsPattern[]] :=
Block @@ Hold[{x, y},
With[{f = Function[, Evaluate[N@f0 /. {x -> #1, y -> #2}], Listable],
stepsize = OptionValue["StepSize"],
step = OptionValue["StepVector"]},
Module[{plot, ptsxy, ptsf, zerodist, pickcrit, contours, split, splitcrit},
plot = ContourPlot[
f @@ ({x, y} + stepsize step) - f @@ ({x, y} - stepsize step) == 0,
{x, x1, x2}, {y, y1, y2},
Evaluate@FilterRules[{opts}, Options[ContourPlot]]];
ptsxy = Cases[plot, GraphicsComplex[pts_, __] :> pts, Infinity];
If[Length@ptsxy > 0,
ptsxy = First@ptsxy;
ptsf = f @@ Transpose @ ptsxy; (* f values at the points *)
zerodist =   (* distances to zero along secant line *)
Abs[ptsf/(f @@
(Transpose @ ptsxy - stepsize step/10) - ptsf /. 0. -> \$MachineEpsilon) *
Norm[stepsize step/10]];
pickcrit = UnitStep[zerodist - OptionValue["Tolerance"]];
contours = Cases[plot, l_Line :> First@l, Infinity];

split = Map[SplitBy[#, pickcrit[[#]] &] &, contours];
splitcrit = Map[pickcrit[[#]] &, split, {-2}];

Graphics[
GraphicsComplex[ptsxy,
Line@Flatten[Pick[split, splitcrit, 0], 1]],
FilterRules[Options[plot], Options[Graphics]]],

Graphics[{},
FilterRules[Options[plot], Options[Graphics]]]
]
]]
];


### OP's first example

I used the default step vector of {0, 1} (because it turned out to work).

zeroSet[Abs@f1[x, y], {x, 0, 4 Pi}, {y, 0, 4 Pi}, PlotPoints -> 80,
MaxRecursion -> 3, PlotRange -> All]


### OP's second example

With[{a0 = 10, a = 14/10},
zeroSet[Abs@f2[a, a0, x, y], {x, 0, 2}, {y, -2 \[Pi]/a0, 2 \[Pi]/a0},
PlotPoints -> 30, MaxRecursion -> 3, PlotRange -> All]
]


Note: The values of f are computed in ContourPlot and again in ptsf. I didn't know how to get the values of f out of ContourPlot in the same order as in the GraphicsComplex. Perhaps there is a way?

Update

As pointed out in the original answer, a limitation of the method is that the function zeroSet does not detect zeros when the curve is (nearly) parallel to the "StepVector". For the example in a comment, one can adjust it.

zeroSet[Abs[Sin[x] Cos[y] Sin[x - y]^2], {x, -3, 3}, {y, -3, 3},
"StepSize" -> 0.001, "Tolerance" -> 0.5, "StepVector" -> {1., -1.},
PlotPoints -> 50, MaxRecursion -> 4]


-
I tried f3[x_, y_] := If[x^2 + y^2 - 4 < 0, 0, x^2 + y^2 - 4]] with zeroSet[f3[x, y], {x, -3, 3}, {y, -3, 3}, PlotPoints -> 80, MaxRecursion -> 3, PlotRange -> All]. I got error messages. –  Hector Sep 26 '13 at 1:20
zeroSet[Abs@f1[x, y], {x, 0, \[Pi]/2}, {y, 0, \[Pi]/2}, PlotPoints -> 80, MaxRecursion -> 3, PlotRange -> All] also chokes but for different reasons. –  Hector Sep 26 '13 at 1:30
zeroSet[Abs[Sin[x] Cos[y] Sin[x - y]^2], {x, -3, 3}, {y, -3, 3}, PlotPoints -> 80, MaxRecursion -> 3, PlotRange -> All] misses some solutions that appear with Module[{lf = Log@ComplexExpand@Abs[Sin[x] Cos[y] Sin[x - y]^2], xp}, xp = Simplify[D[lf, x]^2 + D[lf, y]^2]; ContourPlot[xp, {x, -3, 3}, {y, -3, 3}]]. –  Hector Sep 26 '13 at 3:10

Here is the answer from the Wolfram Technical service, these two methods can produce most of the contours, but still the contours are broken.

Method 1: plot Re[f[k,K0]]+Im[f[k,K0]]==0

a0 = 10.; a = 1.4;
f[k_, K0_] :=
a^2 Sech[(a a0)/2]^2 (-2 I (1 + E^(2 I a0 k)) k +
a (-1 + E^(2 I a0 k)) Tanh[(a a0)/2]) +
2 k (I E^(I a0 k) ((a - k) (a + k) Cos[a0 k] + (a^2 + k^2) Cos[
a0 K0]) + a (-1 + E^(2 I a0 k)) k Tanh[(a a0)/2])

ContourPlot[
Re[f[k, k0]] == -Im[f[k, k0]], {k, 0, 4}, {k0, -2 Pi/a0, 2 Pi/a0},
PlotPoints -> 300, MaxRecursion -> 0,
RegionFunction ->
Function[{k, k0, func}, Exp[Re[f[k, k0]] Im[f[k, k0]]] > 0.99999]]


Method 2: using image processing to combine the contours of the Re and Im parts.

{ContourPlot[Re[f[k, k0]] == 0, {k, 0, 4}, {k0, -2 Pi/a0, 2 Pi/a0},
PlotPoints -> 50, ContourStyle -> Black],
ContourPlot[Im[f[k, k0]] == 0, {k, 0, 4}, {k0, -2 Pi/a0, 2 Pi/a0},
PlotPoints -> 50, ContourStyle -> Black]}


im1 = Image@
ContourPlot[Re[f[k, k0]] == 0, {k, 0, 4}, {k0, -2 Pi/a0, 2 Pi/a0},
ContourStyle -> Black, PlotPoints -> 50, Frame -> None,
im2 = Image@
ContourPlot[Im[f[k, k0]] == 0, {k, 0, 4}, {k0, -2 Pi/a0, 2 Pi/a0},
ContourStyle -> Black, PlotPoints -> 50, Frame -> None,

ImageDimensions[im1]
(*{360, 360}*)
frameticks =
N@{{Transpose[{FindDivisions[{1, 360}, 8],
FindDivisions[{0, 4}, 8]}],
None}, {Transpose[{FindDivisions[{1, 360}, 8],
FindDivisions[{-2 Pi/a0, 2 Pi/a0}, 8]}], None}};
Show[ImageAdd[im1, im2], Frame -> True, FrameTicks -> frameticks]


-

It is a tricky question. CountourPlot[F[x,y]==0,...] finds points where F[x,y]>0 and F[x,y]<0. Then by dichotomy it finds points where F[x,y] is approximately zero. In your case Abs[f[x,y]]-0.001 is almost always positive, so the algorithm fails.

It is nontrivial (at least at the first sight) to find points when both the real and the imaginary part of the complex function is zero. Especially if the function has a quickly oscillating phase.

For this moment I found the following workaround:

It finds points where derivatives of Abs[f[x,y]] have a singularities or are zero.

f[a_, a0_, k_, K0_] :=
a^2 Sech[(a a0)/2]^2 (-2 I (1 + E^(2 I a0 k)) k +
a (-1 + E^(2 I a0 k)) Tanh[(a a0)/2]) +
2 k (I E^(I a0 k) ((a - k) (a + k) Cos[a0 k] + (a^2 + k^2) Cos[
a0 K0]) + a (-1 + E^(2 I a0 k)) k Tanh[(a a0)/2]);
With[{a0 = 10., a = 1.4},
ContourPlot[{Abs[f[a, a0, y, x - 0.01]] == Abs[f[a, a0, y, x + 0.01]],
Abs[f[a, a0, y - 0.01, x]] == Abs[f[a, a0, y + 0.01, x]]},
{x, -\[Pi]/a0, \[Pi]/a0}, {y, 0, 2}]]


In plot only those points are meaningful where both colors match.

-
Do you have any idea why Abs[f[x,y]]==0.001 trick works for (Cos[x] + Cos[y]) Exp[I x y] Exp[x/2] but not my function? –  xslittlegrass Sep 23 '13 at 19:28
@xslittlegrass No, I don't. –  ybeltukov Sep 23 '13 at 20:02