We get 2 solutions:

restart; 
with(Student:-MultivariateCalculus):
A := [0, 0, 0];
B := [c, 0, 0];
S := [x, 0, z];
solve([Distance(S, A) = a, Distance(S, B) = b], [x, z], explicit) assuming a > 0, b > 0, c > 0;

Obviously, the solutions have a geometric meaning if the triangle inequalities  are hold for the parameters  a, b, c . 

Of course, degenerate cases are also possible when the solution is unique ( z=0  and the points  A, B, S  lie on the same line ). It will be if  a=b+c  or c=a+b  or b=a+c :

Sol:=solve([Distance(S, A) = a, Distance(S, B) = b], [x, z], explicit) assuming a > 0, b > 0, c > 0;
simplify(eval(z,Sol[1]), {a=b+c});
simplify(eval(z,Sol[1]), {b=a+c});
simplify(eval(z,Sol[1]), {c=a+b});

the selection from first level operands (as you wanted):

restart;
expr:=x+sin(x)+ln(y)+10+ln(x+y)^2;
{op(expr)}; 
select(t->type(t,{ln(anything),ln(anything)^anything}), %);

                           {ln(x+y)^2, ln(y)}

Edit.

Maybe you should abandon a polar grid and do as in the example below:

restart;
plots:-polarplot(t,t=0..Pi,color=black, thickness=3, coordinateview = [0 .. Pi, 0 .. Pi]);
plot(t,t=0..Pi,color=black, thickness=3, coords=polar, scaling=constrained, gridlines);

Using the data for the solution  C(x,t) , calculate the derivatives numerically at a point  in which  diff(C(x,t),x,x)<>0  and then  d=diff(C(x,t),t)/diff(C(x,t),x,x)

I replaced  I(t)  with  P(t) , removed 2 strange options 0 and reduced the ranges for  P  and  Q  variables for better viewing:

s0 := 3*10^5:
d := 10^(-3):
delta := 10^4:
b := 5*10^6:
lambda := 4.16:

DEtools:-DEplot([diff(P(t), t) = s0 + P(t)*(-d - delta*Q(t)/(b + Q(t))), diff(Q(t), t) = -lambda*Q(t)], [P(t), Q(t)], t = 0 .. 10, P = 0 .. 20,  Q = 0 .. 20, dirfield = 400, arrows = smalltwo, number = 2, [[0, 4, 0.1], [0, 0.2, 4.1], [0, 7, 0.2], [0, 0.2, 7]], color = red, linecolor = blue, numsteps = 100);

expr:=(x^2 - 1)/x^2*y + (x^2 - 1)/x^2;
subs(x^2 - 1 = Z*x^2, expr);

 
Edit.
 

f := a*x+b;
g := unapply(f, a,b,x);
g(a,b,x); 
g(2,3,x); 
g(2,3,t);

I think the following calculation explains this difference:

restart;
ode1:=diff(y(x),x) = exp(x)*sin(x)+y(x)*cot(x);
my_sol:=[sin(x)*(exp(x)+_C1)];
eval(ode1,y(x)=my_sol);

ode2:=diff(y(x),x) = y(x);
my_sol:=[_C1*exp(x)];
eval( ode2,y(x)=my_sol);

That is if  a<>0  then  [a+b] <> a+[b]

Example:

restart;
X:=LinearAlgebra:-RandomMatrix(50,10, generator=-10..10):
Y:=LinearAlgebra:-RandomMatrix(50,10, generator=-10..10):
B1 := plots:-animate(plots:-pointplot, ['X'[round(kk)], 'Y'[round(kk)], symbolsize = 25, symbol = solidcircle, colorscheme = ["valuesplit"]], kk = 1 .. 50, frames=60);

Just apply the  round  command to all the elements of the matrix.

Example:

A:=<1.1,2.5; -3.53,5.7>;
round~(A);

In my opinion, a conceptually simpler code (especially for a beginner) will just go through a double loop on these lists  X  and  Y :

restart;
X := [3, 6, 4, 13]:   Y := [8, 6, 5, 7]:  n:=nops(X):  k:=0:
for i from 1 to n-1 do
for j from i+1 to n do
k:=k+1;  L[k]:=[[[X[i],Y[i]],[X[j],Y[j]]],sqrt((X[i]-X[j])^2+(Y[i]-Y[j])^2)];
od: od:
L:=convert(L, list);
plot(op~(1,L), style= pointline, symbolsize= 20, symbol= solidcircle);

Addition. The code will be 2 times shorter if we use a nested  seq  instead of loops:

restart;
X := [3, 6, 4, 13]:   Y := [8, 6, 5, 7]:  n:=nops(X):  
L:=[seq(seq([[[X[i],Y[i]],[X[j],Y[j]]],sqrt((X[i]-X[j])^2+(Y[i]-Y[j])^2)], j=i+1..n), i=1..n-1)];
plot(op~(1,L), style= pointline, symbolsize= 20, symbol= solidcircle);

Edit.

Try

if  `and`(seq(a[i] <= b[i] , i=1..min(nops(a),nops(b)) ))  then ...  fi;

# Example:
a:=[1,2,3]:
b:=[1,2,4,5]:
if  `and`(seq(a[i] <= b[i] , i=1..min(nops(a),nops(b)) ))  then true else false  fi;

Try  axis = [tickmarks = [7, subticks = 4], thickness = 2]

Edit. At first I did not notice your second plot.  See a workaround below in a corrected file:

Curves_(1)_new.mw

Example:

Download ellipse.mw

The example:

f := exp(2*gamma(t, r)-2*alpha(t, r)-2*beta(t, r));

`*`(exp~([(op@@2)(f)])[]);

This can be applied to any expression  exp(Expr)  if the type of Expr is  `+` .