Get a numerical solution of an equation containing an integral?

I have the following problem:
given r, I have to find the {x,y} which satisfies:

Integrate[1/(E^((2*x^2-4*x*d+2*d^2+2*y^2)/(8*d+4*r^2))*(Sqrt[d]*(2*d+r^2))),{d,0,Infinity}]==1

however, of all the solutions, I am interested to find the one which has the maximum y.
Unfortunately, the integral can be solved only numerically.
Any help?
Bye

If you restrict x, say to be nonnegative, then this is a straightforward optimization problem once the integral is made into a “black box” function that only evaluates for numeric input.

f[x_?NumberQ, y_?NumberQ, r_?NumberQ] := 
 NIntegrate[
   1/(E^((2*x^2 - 4*x*d + 2*d^2 + 2*y^2)/(8*d + 4*r^2))*(Sqrt[
         d]*(2*d + r^2))), {d, 0, Infinity}] - 1

Example:

FindMaximum[{y, f[x, y, 1/2] == 0, x >= 0}, {x, y}]

(* Out[69]= {1.26438464527, {x -> 0.340246814653, y -> 1.26438464527}} *)

Dear Daniel Lichtblau, thank you very much for your double help.
I can not restrict x to be non negative. However, thank to your help, I can now solve the complete problem, that is:

ranger = {0, 0.001, 0.01, 0.1, 1, 10};
rangeRy = {10^-6, 10^-4, 10^-2, 10^0, 10^2, 10^4, 10^6};
f2[x_?NumberQ, y_?NumberQ, z_?NumberQ, r_?NumberQ] := NIntegrate[(1/(E^((x^2 - 2*x*d + d^2 + y^2 + z^2 +(r^2/d)*z^2)/(2*(d + r^2)))*(Sqrt[d]*(d + r^2)))), {d, 0,Infinity}];
solu2 = Table[FindMaximum[{y, Ry/Sqrt[\[Pi]] f2[x, y, 0, r] == 1, y > 0}, {x, y}], {r, ranger}, {Ry, rangeRy}]

As you see, I need to solve it for different r and Ry. Unfortunately, the numerical computation does not seem to converge over the all ranges defined by me, although it easily converge for Ry close to 1, say between 10^-2 - 10^2, and very small values of r, around 0.
You also advised me to reformulate the problem as:

f1[x_, y_, z_, r_] := NIntegrate[(1/(E^((x^2 - 2*x*d + d^2 + y^2 + z^2 + (r^2/d)*z^2)/(2*(d + r^2)))*(Sqrt[d]*(d + r^2)))), {d, 0, Infinity}];
solu1 = Table[NMaximize[{y, Ry/Sqrt[\[Pi]] f1[x, y, 0, r] == 1, y > 0}, {x, y}], {r, ranger}, {Ry, rangeRy}]

but this version has the same problems, even if it seems to converge quicker.

Finally I used the following calculations to check the convergence of the solutions, but often I do not have good results.

MatrixForm[Transpose[Table[rangeRy[[j]]/Sqrt[\[Pi]] f1[solu1[[i, j, 2, 1, 2]], solu1[[i, j, 2, 2, 2]], 0, ranger[[i]]] - 1, {i, 1, Length[ranger]}, {j, 1, Length[rangeRy]}]]]
MatrixForm[Transpose[Table[rangeRy[[j]]/Sqrt[\[Pi]] f2[solu2[[i, j, 2, 1, 2]], solu2[[i, j, 2, 2, 2]], 0, ranger[[i]]] - 1, {i, 1, Length[ranger]}, {j, 1, Length[rangeRy]}]]]

Do I need to increase somehow the accuracy of the FindMaximum or of the NMaximize?

Kind regards
Umberto

scaling analysis 4.nb (23.5 KB)