# Speed up RegionPlot (replacement with ContourPlot?)

I want to plot a region with the following notebook:

Do[{d1[i] = RandomReal[{0, 1}], c1[i] = RandomReal[{0, 2*Pi}]}, {i, 0,
1001}];
Do[a[i] =
Rationalize[Sqrt[d1[i]]*Cos[c1[i]] + I*Sqrt[d1[i]]*Sin[c1[i]],
0], {i, 0, 1001}];
Do[c2[i] = RandomReal[{0, 2*Pi}], {i, 0, 1001}];
Do[b[i] =
Rationalize[
Sqrt[(1 - d1[i])]*Cos[c2[i]] + I*Sqrt[(1 - d1[i])]*Sin[c2[i]],
0], {i, 0, 1001}]

so2[t_, T_, a_, b_, x_,
z_] = (E^(-t - 1/
T) (-b^2 E^(1/T) (2 - 3 E^(t/2) + E^t)^2 x^2 Conjugate[a]^2 +
b (-b E^(1/T)
x ((-8 + 8 E^t - 9 E^(t/2) t + E^(3 t/2) t) (-1 + z) +
E^(1/T) (-8 + 8 E^t - 9 E^(t/2) t) (1 + z)) +
a ((-1 + E^t) (-20 + 24 E^(t/2) + E^t (-4 + t) - 9 t) (-1 +
z)^2 + E^(
2/T) (20 - 24 E^(t/2) + 4 E^t + 9 t) (1 + z)^2 +
2 E^(1/T) (-4 (5 + x^2 - 5 z^2) +
12 E^(t/2) (2 + x^2 - 2 z^2) + 9 t (-1 + z^2) +
6 E^(3 t/2) (-2 + x^2 + 2 z^2) -
E^(2 t) (-2 + x^2 + 2 z^2) +
E^t (6 + 5 t - 13 x^2 - (6 + 5 t) z^2)))) Conjugate[
a] Conjugate[b] +
E^(1/T) Conjugate[
b]^2 (-a^2 (2 - 3 E^(t/2) + E^t)^2 x^2 +
a b x (-(-8 + 8 E^t - 9 E^(t/2) t + E^(3 t/2) t) (-1 + z) -
E^(1/T) (-8 + 8 E^t - 9 E^(t/2) t) (1 + z)) +
b^2 (-4 E^(2/T) (-1 + E^(t/2))^2 (1 + z)^2 -
E^(1/T) (8 - 16 E^(t/2) - 24 E^(3 t/2) + 4 E^(2 t) +
E^t (28 + 9 t)) (-1 + z^2) +
E^t (-1 + z)^2 (-24 - 9 t +
32 Cosh[t/2] + (-8 + t) Cosh[t] + 16 Sinh[t/2] +
t Sinh[t])))))/(2 b (-(-20 + 24 E^(t/2) +
E^t (-4 + t) - 9 t) (-1 + z) +
E^(1/T) (20 - 24 E^(t/2) + 4 E^t + 9 t) (1 + z)) Conjugate[b])

RegionPlot[
Or @@ Flatten[
Table[Re[so2[0.1 + t, 0.4, a[i], b[i], x, z]] < 0, {i, 0, 1000,
1}, {t, 0, 7, 0.1}], 1] && x^2 + z^2 <= 1, {x, -1,
1}, {z, -1, 1}, AxesStyle -> Directive,
LabelStyle -> (FontSize -> 26), FrameLabel -> {x, z},
FrameTicks -> {{{-0.5, 0.5}, None}, {{-0.5, 0.5}, None}}]


The code works, but it is extremely slow. I would like to vary the second parameter of the function so2, so I have to do similar calculation many times.

Is it possible to speed up the calculation?

I have tried to rewrite the code to use the ContourPlot with so2 equals zero, but I have failed. The output is not simply a boundary of RegionPlot, but the
sum of ContourPlots for different i, which does the figure unreadable. What is more the ContourPlot does not seem to work with Table, the code

ContourPlot[
Flatten[Table[
Re[so2[0.1 + t, 0.4, a[i], b[i], x, z]] == 0, {i, 0, 0, 1}, {t, 0,
5, 1}], 1], {x, -1, 1}, {z, -1, 1}, AxesStyle -> Directive,
LabelStyle -> (FontSize -> 26), FrameLabel -> {x, z},
FrameTicks -> {{{-0.5, 0.5}, None}, {{-0.5, 0.5}, None}}]


gives empty output (and I would also like to include x^2 + z^2 = 1 somehow)...

Edit:

I have managed to make ContourPlot with Table work, with

ContourPlot[
Evaluate@Flatten[{Table[
Re[so2[0.1 + t, 0.4, a[i], b[i], x, z]] == 0, {i, 0, 50, 1}, {t,
0, 5, 0.1}], x^2 + z^2 == 1}, 2], {x, -1, 1}, {z, -1, 1},
AxesStyle -> Directive, LabelStyle -> (FontSize -> 26),
FrameLabel -> {x, z},
FrameTicks -> {{{-0.5, 0.5}, None}, {{-0.5, 0.5}, None}}]


which gives the output The code is extremly fast (Timing: 26.1526) in comparision to RegionPlot as given above (with max. i=50 it is still not completed after 10 minutes...).

The problem is: if I make the maximal number i larger, I get an error:

Identity called with 4096 arguments; 1 argument is expected.

What is more from the above picture I would like to extract an RegionPlot as mentioned earlier.