To prevent the automatic change of the argument in the sin, you can use the empty symbol `` .  If any substitutions and other transformations, the expression inside  the extra parentheses retained. To remove these extra parentheses expand command should be applied to the argument:

b:=sin(``(x+y+z-Phi));

c:=eval(%, x=1);

applyop(t->t/2, 1, c);

applyop(expand, 1, c);

First we plot this body:

plots[implicitplot3d](sqrt(x^2+y^2)=5*sin(z)*exp(-(Pi-z)/(1.5)),x=-2.3..2.3,y=-2.3..2.3, z=0..3.7, style=surface, numpoints=1000000,color=khaki, axes=normal, scaling=constrained);

Below the problem is reduced to the calculation of the surface integral over this surface (assuming the surface density equal to 1). Unfortunately Maple cannot symbolically compute the double integral (calculated only one of two successive integrals). The remaining integral is computed numerically:

f:=unapply(5*sin(z)*exp(-2*(Pi-z)/3), z):

F:=unapply(sqrt(f(z)^2-x^2), x,z);

MI_z:=4*int(sqrt(1+diff(F(x,z),x)^2+diff(F(x,z),z)^2)*f(z)^2,[x=0..f(z),z=0..Pi]);

evalf(MI_z);

evalf[20](MI_z);

So the final result   105.1187800

toandhsp, you wrote "How can I add this code into Perpendicularity Criterion?"  You can write your own code (to simplify the equation of the line a*x+b*y+c=0) in the form of a  procedure, and then use this procedure when you want. 

P:=proc(d::equation)

local a, b, k;

a := coeff(lhs(d), x):

b := coeff(lhs(d), y):

k := igcd(a, b):

if a <> 0 then sign(a)*sort(d)/k else sign(b)*sort(d)/k fi;

end proc:

Examples of use:

P(-6*x+2*y+8 = 0);

map(t->[P(t[1]),P(t[2]),t[3]], [[-4*x-3*y-44 = 0, -4*x+3*y-14 = 0, [x = -29/4, y = -5]], [-4*x-3*y-44 = 0, -3*x-4*y-48 = 0, [x = -32/7, y = -60/7]], [-4*x-3*y-44 = 0, 3*x+4*y-2 = 0, [x = -26, y = 20]]]);

                                                                    3 x - y - 4 = 0

[[4*x+3*y+44 = 0, 4*x-3*y+14 = 0, [x = -29/4, y = -5]], [4*x+3*y+44 = 0, 3*x+4*y+48 = 0, [x = -32/7, y = -60/7]], [4*x+3*y+44 = 0, 3*x+4*y-2 = 0, [x = -26, y = 20]]]

psi:=u(x,y,t)+I*v(x,y,t):

evalc(diff(psi,t)*I+diff(psi,x,x)=epsilon*diff(evalc(abs(psi)^2),y,y)*psi);

evalc~([(Re,Im)~(%)]);  # This is your desired system

Edited.

You can display it as follows:

aa1:=``(f(x+Dx,y+Dy)-f(x,y+Dy)) - ``(f(x,y)-f(x,y+Dy));

Obviously that in real domain the system of three equations  {pd1=0, pd2=0, pd3=0}  is equivalent to the single equation  pd1^2+ pd2^2+ pd3^2=0

restart;

AllGraphs:=proc(n::posint)

local V, S;

uses combinat, GraphTheory;

V:={$1..n};

choose(V,2);

S:=powerset(%);

[seq(Graph(n,s), s=S)];

end proc:

Examples of use:

AllGraphs(3);

plots[display](Matrix(2,4, GraphTheory[DrawGraph]~(%)));

or

AllGraphs(4);

plots[display](Matrix(8,8, GraphTheory[DrawGraph]~(%)));

F:=(d^4-2)*C+(7*d^3-3*d)*C^2-(10*d^4-4*d)*L^2+(d-d^2)*L^3+(R+z^2)*x1+(10*d^3-4*d)*L:

F:=coeff(F,C)*C+coeff(F,L)*L+coeff(F,x1)*x1;

                                          F := (d^4-2)*C+(10*d^3-4*d)*L+(z^2+R)*x1

or more automatically

restart;

F:=(d^4-2)*C+(7*d^3-3*d)*C^2-(10*d^4-4*d)*L^2+(d-d^2)*L^3+(R+z^2)*x1+(10*d^3-4*d)*L:

V:=<C,L,x1>:

F:=add(coeff(F,v)*v, v=V);

Because the function has circular symmetry, it is natural to immediately use the polar coordinates:

int(int(Heaviside(1-r^2)*r, r=0..1), phi=0..2*Pi);

                                      Pi

This calculation can be easily done by hand (without Maple).

Student[Basics][ExpandSteps](factor(x^4+4));

You have to read from the bottom up.

display(fig1, plot(eval(Model, MinSol[2]), x = 0 .. 1.5, color = red, thickness = 0));

You can generate a random function of a certain class of functions depending on a finite number of parameters.

Example: a random function of the class  a*sin(b*x+c)  in which  [a,b,c]  is a random list of floats with the uniform distribution in the range  0..10 :

RandomTools[Generate](list(float(range=0...10, method=uniform), 3));

eval(a*sin(b*x+c), [a,b,c]=~%);

Should be  D(x)(0)  instead of  D(x) (0)

A space between  D(x)  and  (0)  Maple interprets as a multiplication.

For any real numbers b>0, c>0  the expression sqrt(sqrt(9)*sqrt((1+(b+1)^2*c^2+((10/3)*b-2)*c)*(1+(b+1)^2*c^2+(2*b-2)*c))+3+(3*b^2+6*b+3)*c^2+(8*b-6)*c)  will be strictly positive . Here is the simple solution (in fact by hand):

u:=sqrt(sqrt(9)*sqrt((1+(b+1)^2*c^2+((10/3)*b-2)*c)*(1+(b+1)^2*c^2+(2*b-2)*c))+3+(3*b^2+6*b+3)*c^2+(8*b-6)*c);  # The original expression

v:=normal(u^2);  # The simplified expression under the square root

v1:=3*b^2*c^2+6*b*c^2+3^(1/2)*((3*b^2*c^2+6*b*c^2+10*b*c+3*(c-1)^2)*(b^2*c^2+2*b*c^2+2*b*c+(c-1)^2))^(1/2)+8*c*b+3*(c-1)^2;  # In  v  we substitute c^2-2*c+1=(c-1)^2  # All terms are strictly positive if  c<>1, c>0, b>0

eval(v1, c=1);  # for c=1

If the parameters  beta, gamma, R4 and C2  are not necessarily positive, vv's solution does not work. Assume  command helps in this case. Also I wrote  local gamma;  because  gamma is protected constant (Euler's constant).

restart;

local gamma;

assume(gamma>0,beta<0,R4>0,C2<0):

f:=gamma*sqrt(4)*sqrt(omega^2*C2^2*R4^2/(C2^4*R4^4*beta^2*gamma^2*omega^4+C2^2*(1+gamma^2*(beta+1)^2-2*gamma)*R4^2*omega^2+1)):

diff(f,omega):

solve(%=0, omega);

select(t->is(t>0),[%])[];