sqrt(-2*cos(alpha)*cos(alpha+d)+2-2*sin(alpha+d)*sin(alpha))/d;
combine(%);
#MultiSeries:-
limit(%, d = 0, right);

                                         1/2
                          (-2 cos(d) + 2)
                          ------------------
                                  d


                                  1

One can clean up the task (and use g instead of gamma = Euler gamma)

subs(gamma=g, D=d*sin(theta), integrand); # theta <> k*Pi
A:=collect(%,[g, cos(psi), sin(theta)]);
A0:=(1/(-sin(phi)*r+d)-1/(sin(phi)*r+d))*r^2*sin(phi);

'A=2*cos(theta)/sin(theta)*g*cos(psi) * A0'; is(%);

                                 true

eval(A0, d=r*a);
A1:=collect(%, r);

                 /      1              1      \
           A1 := |------------- - ------------| sin(phi) r
                 \-sin(phi) + a   sin(phi) + a/


Integrating w.r.t. r gives aboveTrig * L^2/2 and it remains to find

L^2/2 * 2*Int((1/(-sin(phi)+a)-1/(sin(phi)+a))*sin(phi),phi = 0 .. 1/2*Pi)

Simplify + value gives it, simplify further if desired.

I do not discuss the cases when the above substitutions are invalid
due to zero values (nor checked when Maple uses 'generic' assumptions).

It works as well for the other way round, i.e. first w.r.t. phi.

expand(eq2(x,t)*cosh(1.8751*x/l));
Int(%, x = 0 .. l);
Expand(%);
value(%);

40.0030457235592*(diff(diff(q[1](t), t), t))+201.223915538476*q[1](t)+(30022.0388740059+0.417825675567801e-13*I)*q[1](t)*q[4](t)^2-135846.229059406

f(x)=(2.25*f(x-1)-0.5*f(x-2));
convert(%, rational);
rsolve({%, f(1)=1/3, f(2)=1/12}, f(t));

PS: is it homework?

100*sqrt(a) = .7559289460*sqrt(a)-2.000000000;
isolate(%, sqrt(a));

                       1/2
                      a    = -0.0201523373513343

Up to a 'scaling it is the same as sqrt(x) = - 1. Now take squares on both sides,
sqrt(x)*sqrt(x) = (- 1)*(-1) and you see why.

I agree with Markiyan, it is "illposed" --> Wikipedia. 

You have a polynomial of total degree ~ 3000 and it not likely that
usual precision is enough (at least I would not trust Matlab at all,
except if you are sure that evaluation is stabilized).

Using "ldegree" for your variables suggest to divide by factor1 := 
  phi[1]^45*phi[2]^49*phi[3]^64*phi[4]^62*
  phi[5]^44*phi[6]^36*phi[7]^51*phi[8]^48

  for sing in Vars do
    ldeg, hdeg = ldegree(tarfun, sing), degree(tarfun, sing);
    print(sing, %);
  end do: sing:='sing':

Now consider T1 := tarfun/factor1

NB: do not try to expand to recognize it as polynomial ...

Then you still have degrees of ~ 150 in the most variables while for
eta[p2], mu[p] you have degree = 1366 (try to feed Markiyan's result
for that [hm ... why is eta[p2] negative? It should be in 0 .. 1 ?]).

For the coefficients perhaps it may help to scale T1 * 10^(-472).

NB: I think that these numerical problems remain even if you take
log of the (modified) function.

Then you may 'review' you constraints to reduce its occurencies, if
that is possible - they roughly describe 'linear diagonal conditions'
in the hypercube [0 .. 1]^18.

There are approaches where those constarints introduce new variables
to be solved for (thus they increase the problem).

Finally: I do not believe that you can 'solve' it by brute force, may
be you need to study the geometry of the actual task.

The integrand ln(sec(x)+tan(x)) = ln((sin(x)+1)/cos(x)) has period = 2*Pi,
so one period is enough. 

We start in x = -Pi/2 and Change(%, x=xi+1/2*Pi, xi): subs(xi=x, %);
makes it live over -Pi .. Pi, Int(ln(-(1+cos(x))/sin(x)),x=-Pi .. t )
then is the new task.

Note that the integral exist: using series the problem for that is
ln(1/x) to be integrated, which can be done.

  Int(ln(-(1+cos(x))/sin(x)),x=-Pi .. t, epsilon=1e-4 );
  plot([Re(%),Im(%)], t= -Pi .. Pi, color=[red,blue]);

shows a symmetry for the real part and a linear imaginary part.

This can be understood by plotting -(1+cos(x))/sin(x), its symmetry
and using ln(-y) = ln(y) + Pi*I (for 0 < y). The last summand gives
the linear imaginary part.

Now restrict to -Pi <= t <= 0 and evaluate the integral symbolically.

Plotting is a mess beyond - but fine for that range and indicates that
it is continous, imaginary = constant (I have not examinated further).

We only need the real part (due to the range) and the constant.

After tedious trials (needing Heaviside to kill derivative signum)
finally for that range I came up with 

sol := 
  -Im(dilog(1/(-cos(t)+1)*sin(t)*I)+dilog(1/(-cos(t)+1)*sin(t)*I+1)) +
    1/2*ln(1/(-cos(t)+1)*(cos(t)+1))*arctan(-1/(-cos(t)+1)*sin(t)) +
    1/4*Pi*ln(2)+1/4*Pi*ln(1/(-cos(t)+1));

sol = Int(ln(-(1+cos(x))/sin(x)), x = -Pi .. t);

For the full range -Pi <= t <= +Pi then one uses

piecewise(-Pi < t and t <= 0,   sol, 
            0 < t and t < +Pi, -eval(sol, t=-t) + Pi*I*t);

For the 'boundary' in t = +-Pi or t=0 one evaluates by taking the limit.

My worksheet is currently to ugly to post it.

Quite often --> hypergeom --> Standardfunctions works:

  Int(ln(1+x)/(1+x^2),x=0..1);  value(%);
  convert(%, hypergeom);
  convert(%, StandardFunctions);
  simplify(%);

                             1/8 Pi ln(2)

sqrt(...) is complex valued on negative inputs. You can - for example - try:

int((g2(t)-g1(t))^2, t=0.. 1);
J:=simplify(%);

Optimization[Minimize](sqrt( abs(J) ));

Edit: you can do it without abs or sqrt

Optimization[Minimize]( J );
map(fnormal, %);
map(identify, %);
lprint(%);
sol:=%[2];

    sol := [c[1,0] = 0, c[1,1] = 1/6, c[1,2] = 1/3, c[2,0] = 1/3, c[2,1] = 1/2, c[2,2] = 2/3, 
              c[3,0] = 2/3, c[3,1] = 5/6, c[3,2] = 1]

eval(J, sol);
                                  0

[g2(t),g1(t)]: eval(%, sol): simplify(%);

                            { 0        t < 0
                            {
                        [t, { t        t < 1 ]
                            {
                            { 0        1 <= t

plot(%, t=0..1, color=[blue,red],thickness=[1,3],discont); # a straight line

[g2(t) - g1(t)]: eval(%, sol): simplify(%);

                          { t        t < 0
                          {
                         [{ 0        t < 1 ]
                          {
                          { t        1 <= t

plot(%, t=0..1, color=[blue,red],thickness=[1,3],discont); # zero, except in t=1

That antiderivative is not continous, you have to split (where - plot and 'guess') and taking the real part should do

Edited: look at your integrand, it has imaginary parts

  ln((1+sin(x))/cos(x));
  plot([Re(%),Im(%)],x=0 .. 10);

I get 0, 3.45314652776392*I, +-2.11652746155732 -.304189592310899*I, but
there may be more as the task is numerically very ugly.

Sketch: 

Taking the series in 0 one sees to scale by const := -.128305934050000e27
that omega = 0 is fine and to consider

  Frequency_Equation/const:  
  map(combine, %, power): simplify(%):collect(%, omega): simplify(%, size):
  #convert(%, rational):
  E:=%;

Now have a look at it it around the origin

  theRange:=  x=-2*Pi .. 2*Pi , y=-2*Pi .. 2*Pi;
  eval(E, omega=x+I*y);
  plot3d(abs(%),  theRange, view= 0 .. 3, axes=boxed, orientation=[-40,60]);

From that look at theRange:=  x=-0.1 .. 2.5 , y=0 .. -1;
indicating to take

  fsolve(E, omega= +2 - 0.3*I, complex);
  fsolve(E, omega= -2 - 0.3*I, complex);

giving the last 2 solutions.

Now the spurious paek in the first plot also might be a solution, try
theRange:=  x=-0.1 .. 0.1 , y=3 .. 4; and 

  fsolve(E, omega= 0 +3.5*I, complex);

gives the remaining one.

Note that it has to be taken with care due to ugly numerical situation.

Here is a way with usual precision (having Digits = 15, but 10 would do).

  f := (a, b, k) -> exp(-x)^k*(Sum(x^m/(m!), m = a .. b))^k; # exp outside

  h1 := f(0,0,82)*f(1,3,49)*f(4,6,47)*f(7,10,47)*f(11,15,57)*f(16,20,40)
        *f(21,25,38)*f(26,35,52);
  h0:=(1-f(0,35,1));
  h2:=h0^91;

Then h = h1*h2. For 0 < x it is clear that h1 has positive factors and
for h0 that follows from Sum(x^m/m!,m = 0 .. 35) <= exp(x) = Sum to oo.

Now apply log to show the problem and solve it (log preserves maximum):

  ln(h1): simplify(%) assuming 0<x;
  LNH:= % + 91*ln(h0);       # helping Maple
  plot(%/ln(10), x=6 .. 20); # to show decimal magnitude --> 10^(1560), acer


  diff(LNH, x): value(%):
  fsolve(%, x=15);
                           15.6881165580849

Int(5*exp(t)*exp(-2*j*pi*f*t), t = -infinity .. 0)+
Int(5*exp(-t)*exp(-2*j*pi*f*t), t = 0 .. infinity);

subs(j=I, pi=Pi, %); value(%); simplify(%);

                                  10
                             ------------
                                 2  2
                             4 Pi  f  + 1

Lower case sum enforces to find symbolic solutions for each call

You may also consider something like

F:= (t,M) -> Sum(sin(n*t), n=0..M);
plot( F(t,11), t=0..2*Pi);

PS: please post input as text instead of pictures (or at additionally provide it), nobody wants to type in your task

I guess: x^(x^x) + floor( x^(x+1)*ln(x)^4 ), x=10 (floor -> 2811012357389)
But have no sound proof (and what I have seen was to tedious for me,
the keywords are "tetration" and "iterated exponentiation").

Roughly: write x^(x^(x^(x^r))) as series in r, substract x^(x^x) and now
set r=x^(-x^x). 'Simplify', then abbreviate by x^(x^x) = xi. Feed x=10.

This means to divide by xi = extremely large, so floor only lives on the
term x^(x+1)*ln(x)^4, giving the asserted: rest is very small positive,
smaller than frac of that, ~ 0.44

For numerics first apply log, then feed x=10.0 (as float), take exp and
apply floor (and frac).

Using series in that way is not sound, but I think it could be justified.