You mean to construct a right circular cone given radius and height?  It's very simple:

R:=1:  H:=2:

A:=plot3d([r*cos(phi), r*sin(phi), H-H/R*r], r=0..R, phi=0..2*Pi):  # The lateral surface

B:=plot3d([r*cos(phi), r*sin(phi), 0], r=0..R, phi=0..2*Pi):  # The bottom

plots[display](A,B, axes=normal, scaling=constrained);

Very interesting integral! I have only one question to Markiyan. Where did you find it?

The basic idea of ​​the solution - two-stage estimates of the integrand at every interval  x=Pi*k/2 .. Pi*(k+1)/2 , k=1..infinity . Of course, all of the estimates made ​​by hand. Maple was used only for presentation of results, visualizations and calculation  of the integral on every interval.

We have the following obvious estimates for every interval

The latest estimate holds for odd half-intervals  x=Pi*k/2 .. Pi*k/2+Pi/4, as  

sin(x)^2 >= cos(x)^2, 2*sin(x)^2 >= sin(x)^2+cos(x)^2, 2*sin(x)^2 <= -2*x+(1/2)*Pi+Pi*k+1

For even half-intervals all similar, but  sin  should be replaced by cos . 

Visualization of the estimates  2*sin(x)^2 <= -2*x+(1/2)*Pi+Pi*k+1  and  2*cos(x)^2 <= -2*x+(1/2)*Pi+Pi*k+1 for  x=Pi/2..3*Pi/2 :

plot([2*sin(x)^2, 2*cos(x)^2, -2*x+3*Pi/2+1,-2*x+5*Pi/2+1], x=Pi/2..3*Pi/2, 0..Pi/2+1, thickness=[1,1,2,2], color=[red,blue,red,blue]);

Visualization of the chain of inequalities for 3 intervals: 

f:=x->1/(x^(2*sin(x)^2)+x^(2*cos(x)^2)):

g:=x->piecewise(x>Pi/2 and x<=Pi,1/(Pi^(2*sin(x)^2)+Pi^(2*cos(x)^2)), x>=Pi and x<=3*Pi/2, 1/((3*Pi/2)^(2*sin(x)^2)+(3*Pi/2)^(2*cos(x)^2)), x>=3*Pi/2 and x<=5*Pi/2, 1/((5*Pi/2)^(2*sin(x)^2)+(5*Pi/2)^(2*cos(x)^2))):

h:=x->piecewise(x>Pi/2 and x<=3*Pi/4, 1/2/Pi^(-2*x+3*Pi/2+1), x>Pi and x<=5*Pi/4, 1/2/(3*Pi/2)^(-2*x+5*Pi/2+1),  x>3*Pi/2 and x<=7*Pi/4, 1/2/(5*Pi/2)^(-2*x+7*Pi/2+1)):

plot([h, g, f], Pi/2..2*Pi, color=[red,blue,green], filled);

Next we compute integrals of the function   h(x)  (red shapes) on every  half-interval and estimate them:

int(1/(2*(Pi*(k+1)/2)^(-2*x+Pi/2+Pi*k+1)), x=Pi*k/2..Pi*k/2+Pi/4);

The series 

 diverges.

for k from 0 to 629 do

if msolve(119^x = k, 630)<>NULL then print(['k' = k, msolve(119^x = k, 630)])

end if; end do;

We are satisfied only  k = 91, because it is obvious that the number  30^(8^33)  is divisible by 6

In my opinion both calculations are correct. Regarding the former, see help on  LinearAlgebra[Add]  command. In the second calculation  a  is a unknown object, it can be either scalar or matrix, so Maple leaves this expression unevaluated.

for  first  fraction:

restart;

a:=ifactors(numer(convert(0.999987406876435, fraction)));

b:=ifactors(denom(convert(0.999987406876435, fraction)));

n:=nops(a[2]): m:=nops(b[2]):

frac1:=convert([seq((x||i)^a[2,i,2], i=1..n)], `*`)/convert([seq((x||i)^b[2,i-n,2], i=n+1..m+n)], `*`);

assign(seq(x||i=a[2,i,1], i=1..n), seq(x||i=b[2,i-n,1], i=n+1..m+n));

seq(x||i, i=1..m+n):

seq(cat('x', i)=x||i, i=1..m+n);

a:=ifactor(op(1, convert(0.999987406876435, fraction)));

b:=ifactor(op(2, convert(0.999987406876435, fraction)));

c:=ifactor(op(1, convert(0.999919848203811, fraction)));

d:=ifactor(op(2, convert(0.999919848203811, fraction))); 

x1:=2:

x2:=3:

x3:=7:

x4:=17:

x5:=173:

x6:=709:

x7:=5347:

x8:=18713:

L:=[seq([x||i, x[i]], i=1..8)]:

for k in [a, b, c, d] do

subs({seq(L[i,1]=cat('x', op(L[i,2])), i=1..8)}, k);

od;

convert / phaseamp  appeared only in recent versions of Maple. If you have an older version, you can use  applyrule  command:

Rule:=a::realcons*sin(x)+b::realcons*cos(x)=sqrt(a^2+b^2)*sin(x+arctan(b,a)):

applyrule(Rule, sin(x)+cos(x));

To plot a graph of  phi  against time is  easier than making of animation. This can be done directly by the command  plots:-odeplot  using the procedure  integ :

ode := {diff(phi(t), t, t) = -phi(t), phi(0) = (1/12)*Pi, (D(phi))(0) = 0};

integ := dsolve(ode, numeric);

plots:-odeplot(integ, [t, phi(t)], 0..15);

We assume that  A - is a known matrix, for example,  A = <<1 | -2 | 3> , <4 | 1 | 2> , <5 | 1 | 2>> ,  and  B - is a subspace of  R^3, for example, all vectors collinear to  <1, 2, 3> . First, we find the general form of the matrices  Y  for which  NullSpace(Y) = B :

restart;

Y:=<<a[1] | b[1] | c[1]> , <a[2] | b[2] | c[2]> , <a[3] | b[3] | c[3]>>:

solve({op(Equate(Y.<1, 2, 3>, <0, 0, 0>))}, {a[1], a[2], a[3]});  assign(%);

Y;

Columns  <b[1], b[2], b[3]>  and  <c[1], c[2], c[3]>  should not be proportional.

Next, we find the matrix  X :

X:=Matrix(3, symbol=x);

A:=<<1 | -2 | 3> , <4 | 1 | 2> , <5 | 1 | 2>>;

Sys:=op(Equate(A.X+X.A, Y));

solve({Sys}, {seq(seq(x[i,j], j=1..3), i=1..3)});

We see that there are infinitely many solutions, depending on the 6 parameters.

You can use  RealRange command . For example:

is(2 in RealRange(0, Open(2)));

is(1 in RealRange(0, Open(2)));

is(0 in RealRange(0, Open(2)));

                              false

                              true

                              true

When solving equations, you can just use the corresponding inequalities. Because in your equation 3 unknowns, two of them must be specified. For example:

restart;

y:=2:  z:=0.3:

solve({y=1.048 + 1.02*x + 6.118*(z-4.041*x^2) + 16.22*(z^2) +6.241* (x*z),

x>=0.001, x<=0.7}, x);

                                                         {x =0 .3718638796}

a := x->piecewise(x^2 <= 3, x^2, x^1.5):

a(1);  a(2);

                  1

        2.828427125

Addition: one more way

a := x->`if`(x^2 <= 3, x^2, x^1.5):

convert(exp(x), Sum, dummy=n);

This is quite easy. For example:

V:=[V[k] $ k=1..n]: N:=[$ 1..n]:

plot(N, V, thickness=2);

Even for a specified  b the  solution is expressed only in terms of integrals:

for b from -5 to 5 do

P[b](t):=t*ln(t)^b:

dsolve(diff(Q(t), t)^2 = diff(P[b](t), t), Q(t)):

od;