/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 16.04.0 ] */

/* [wxMaxima: input   start ] */
kill(all)$fpprec:32$fpprintprec:6$ratprint:false$load(draw)$set_draw_defaults(terminal=wxt,grid=true)$
N:9$  wp:1$  LP:1$
a[k]:=(2*N-k)!/(2^(N-k)*k!*(N-k)!)$ A[k]:=2*(N-k)/((2*N-k)*(k+1))$
Pi(s):=sum(a[k]*s^(k),k,0,N)$ Pi(s);   Q(s):=sum(a[k]/a[0]^((N-k)/N)*s^k,k,0,N)$ Q(s),numer;
bfallroots(Pi(s))$ S:map(rhs,%)$  ss:sort(S)$
re(x):=realpart(S[x])$ im(x):=imagpart(S[x])$
B(s):=a[0]/Pi(s)$
for i:1 thru N step 2 do print('re[(i+1)/2],"=",
 if oddp(N+1) and im(i)=0 then [float(re(i)/wp),float(re(i+1)/wp)] else float(re(i)/wp)," , ",
 'im[(i+1)/2],"=",float(im(i)/wp) )$
draw2d(point_type=5,point_size=0.25,/*proportional_axes=xy,*/points_joined=true,yrange=[0,N-1],
    points(realpart(ss),imagpart(ss)),
    color=red,nticks=100,parametric(re(1)*(1+(cosh(t*1.5)-1)/(1-cosh(1.5)))^1,(N+0.5)*t,t,-1,1) )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
sx:sort(abs(imagpart(S)))$ d:makelist(sx[k+1]-sx[k],k,1,length(S)-1,2)$
wxplot2d([discrete,d]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
105/2.11394**4;15/1.75569**3;3/1.36165**[1,2];%phi,numer;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
H(s):=a[0]/(prod(s^2+2*abs(re(n))*s+re(n)^2+im(n)^2,n,1,ceiling(N/2))*(if oddp(N) then s+abs(re((N+1)/2)) else 1))$ H(s);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
for i:1 thru N-1 do print([a[i+1]/a[i],A[i]]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
float(945**.2);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
HP(s):=a[0]/sum( a[k]/s^k,k,0,N )$ ratsimp(HP(s));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
kill(all)$ fpprinprec:6$ ratprint:false$
for n:1b0 thru 32b0 do [ Den(s):=sum((2*n-k)!/(2^(n-k)*k!*(n-k)!)*s^(k),k,0,n),r:allroots(Den(s)=0),
 real(k):=realpart(rhs(r[k])),imag(k):=imagpart(rhs(r[k])), /* for i:1 thru n step 2 do */
  for j:1 thru n step 2 do ( ([n,(j-1)/2+1],if j<=n then -[real(j),imag(j)] else 0 ), n,
   sprint(j,-real(j)),if j=n then newline() ) ]$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
  [real(1),real(2),real(3),real(4),real(5),real(6),real(7),real(8),real(9),real(10),
   real(11),real(12),real(13),real(14),real(15),real(16),real(17),real(18),real(19),real(20),
   real(21),real(22),real(23),real(24),real(25)]
   [wxMaxima: comment end   ] */

/* [wxMaxima: hide output   ] */
/* [wxMaxima: input   start ] */
for n:1 thru N do [ Den(s):=sum((2*n-k)!/(2^(n-k)*k!*(n-k)!)*s^(k),k,0,n),r:allroots(Den(s)=0), j:2,
 print(-if oddp(n+1) and imagpart(rhs(r[j]))=0 then realpart([rhs(r[j]),rhs(r[j+1])]) else realpart(rhs(r[j])) ) ]$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
X(s):=1.0*s^25+325.0*s^24+52650.0*s^23+5651100.*s^22+4.506752*10^8*s^21+2.83925*10^10*s^20+1.46695*10^12*s^19+6.37074*10^13*s^18+2.36514*10^15*s^17
+7.59472*10^16*s^16+2.12652*10^18*s^15+5.21965*10^19*s^14+1.12657*10^21*s^13+2.14049*10^22*s^12+3.57768*10^23*s^11+5.24726*10^24*s^10+6.72305*10^25*s^9
+7.47445*10^26*s^8+7.14225*10^27*s^7+5.78898*10^28*s^6+3.90756*10^29*s^5+2.13986*10^30*s^4+9.14302*10^30*s^3+2.86216*10^31*s^2+5.84358*10^31*s+
5.84358*10^31$
allroots(X(s)=0);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
par(s):=sum(a[2*k]*(%i*s)^(2*k),k,0,N/2)$ imp(s):=sum(a[2*k-1]*(%i*s)^(2*k-1),k,1,ceiling(N/2))$ [par(s),imp(s)];
RE:sum( coeff(par(s),s^k),k,0,N )+a[0]$  IM:sum( coeff(imp(s),s^(2*k+1)),k,0,ceiling(N/2) )$ [RE,IM];
X:a[0]/sqrt(RE^2+IM^2)$ Att:-20*log(X)/log(10)$ float([X,1/X,1/2**.5*X]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
poli(s):=sqrt((par(s))^2+(imp(s))^2)$ float(rectform(-20*log(poli(8.75))/log(10)));
/* wxplot2d([(a[0]/poli(%i*x/2/%pi))],[x,.1,10],[logx],[logy],[gnuplot_preamble,"set grid"],[grid,5,5])$ */
float(poli(8.75)^2);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
re(x):=abs(rhs(realpart(S(k)[x])))$  im(x):=abs(rhs(imagpart(S(k)[x])))$
for i:2 thru N step 2 do print( sigma[i/2],"=",float(re(i))," , ",omega[i/2],"=",float(im(i)) )$
if oddp(N) then print( sigma[(N+1)/2],"=",float(re(N) ) )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
rehp(x):=re(x)/(re(x)^2+im(x)^2)$   imhp(x):=im(x)/(re(x)^2+im(x)^2)$
for i:2 thru N step 2 do print( sigma[i/2],"=",float(rehp(i))," , ",omega[i/2],"=",float(imhp(i)) )$
if oddp(N) then print( sigma[(N+1)/2],"=",float(rehp(N) ) )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
a_2:0$  o(x):=(2*(N-(1-a_2)*N+signum(a_2-0.5)*x)/(2*N-(1-a_2)*N+signum(a_2-0.5)*x)/
             ((1-a_2)*N+signum(a_2-0.5)*x+1))^signum(a_2-0.5)$
orig(x):=( (2*N-((N-1)*a_2-signum(a_2-0.5)*x))*(((N-1)*a_2-signum(a_2-0.5)*x)+1)/2/
             (N-((N-1)*a_2-signum(a_2-0.5)*x)) )^(signum(a_2-0.5));
qq(x):= ((2*N-(a_2*N+a_2-1+((-1)**a_2)*x))*((a_2*N+a_2-1+((-1)**a_2)*x)+1)/2/(N-(a_2*N+a_2-1+((-1)**a_2)*x)))**((-1)**a_2)$
q(x):= ((2*N*(1-a_2)+2*a_2*N-x)*(a_2*N-signum(a_2-0.5)*x+1)/2/(N*(1-a_2)+N*a_2-x))**((-1)**0)$
for i:0 thru N do print(float( [prod(orig(j),j,0,i-1)/105,a[N-i]/105] ))$ float(orig(N-1));
for i:0 thru N step 1 do ( float( 1*prod(o(j),j,1,i) ) )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
Ph(s):=s^N/sum( a[k]/a[0]*s^(N-k), k, 0, N )$ float(Ph(s));
P2(s):=s^(N*a_2)/sum( prod(orig(j)*s,j,0,k-1)/a[0]^a_2, k, 0, N )$ fullratsimp(P2(s),s);
Bx(s):=B(1/s)$ fullratsimp(Bx(s));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
X(k):=(2*N-k)*(k+1)/(2*(N-k)); b(x):=prod( X(N-k), k, 1, x )$
for k:0 thru N do print( float(b(k))," , ", float(X(k)) )$
A(s):=sum( (%i)^k*a(k)*s^k,k,0,N )$ A(s);
Ar(s):=sum( %i^(2*k)*a(2*k)*s^(2*k),k,0,N/2 )$ Ar(s);
Ai(s):=sum( %i^(2*k+1)*a(2*k+1)*s^(2*k+1),k,0,(N-1)/2 )$ Ai(s);
H(s):=a(0)/sqrt( realpart(A(s))^2+imagpart(A(s))^2 )$
H2(s):=a(0)/sqrt( Ar(s)^2+Ai(s)^2 )$
float( H(1) );
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
for w:N/3 step .0001 do  if float(H(w))<=float(10^(-3/20)) then return(w);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
plot2d([H(x)],[x,.01,100],[logx],[gnuplot_preamble, "set grid"]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
(B(0.3556)); B(.35066); float(B(1));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
1/.359381;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
B[0](s):=1$  B[1](s):=s+1$
for i:2 thru N do B[i](s):=(2*i-1)*B[i-1](s)+s^2*B[i-2](s)$
factor(B[N](s));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
for k:0 thru N do print([[(2*N-k)!,2^(N-k),k!,(N-k)!],[a(k),k]])$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
H(s):=Pi(N)$ ss:%i*omega$ zq:-3$
r(s):=realpart(H(s))$ i(s):=imagpart(H(s))$
Phi(s):=atan(i(s)/r(s))$ A(s):=20*log(a(0)/sqrt(r(s)^2+i(s)^2))/log(10)$
for omega:0.1 thru 1.6 step 0.1 do ('Phi,(s)[omega],"=",float(Phi(%i*omega)))$
inc:0.0001$ guess:if N=2 then 1.35 else if N=3 then 1.74 else if N=4 then 2.1 else
if N=5 then 2.41 else if N=6 then 2.69 else if N=7 then 2.94 else if N=8 then 3.16
else if N=9 then 3.38 else if N=10 then 3.57 else if N=11 then 3.74 else 3.94$
for sc:3.58 step inc while A(%i*sc)>=zq
 do print(if not A(%i*(sc+inc))>=zq then
 (scale:sc, print('A(s),"=",float(A(%i*sc))," ; ",'𝜔,"=",float(sc),
 " ; ",'A(s+inc),"=",float(A(%i*(sc+inc)))," ; ",'𝜔,"=",float(sc+inc))) ) $
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
kill(sc);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
inc2:0.1$ scale:3$
for n:1 thru N do
 for a:0.1 step 0.1 thru -dB do
  for sc2:scale step inc do
   while A(%i*sc2)<=a do
    if not A(%i*(sc2+inc))>a then
     print(float(sc2),float(a)),scale:sc2,
      for sc2:scale step inc2 do
       while A(%i*sc2)<=a do
        if not A(%i*(sc2+inc2))>a then
         (print('A(s),"=",float(A(%i*sc2))," ; ",omega,"=",float(sc2)),sc2:float(sc2),float(a) );
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
re(k):=realpart(S(k))$ im(k):=imagpart(S(k))$
rh(k):=re(k)/(re(k)^2+im(k)^2)$ ih(k):=im(k)/(re(k)^2+im(k)^2)$
for k:1 thru N/2 step 2 do ('𝕽,'e[k],"=",float(re(k))," ; ",'𝕴,m[k],"=",float(im(k)))$
if LP=1 then print(float(re(k))," ; ",float(im(k))) else print(float(rh(k))," ; ",float(ih(k)))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
if LP=1 then print('ω[0]^2,"=",float(re(k)^2+im(k)^2)," ; ",'2,"*",'𝕽,'e,"=",float(2*re(k))," ; ")
else print('ω[0]^2,"=",float(rh(k)^2+ih(k)^2)," ; ",'2,"*",'𝕽,'e,"=",float(2*rh(k))," ; ")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
fs(x,y):=a*exp(-0.5*(((log(x)-b)/c)^2 + ((log(y)-d)/f)^2)) + Offset$
a :  8.2817504399961680E+06$
b :  5.4666920855164321E+01$
c :  1.0465433023210867E+01$
d :  8.1370007073912198E+00$
f :  3.3040744355996408E+00$
Offset :  2.0569866296476450E-02$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
fs(x,y):=a*x^b*y^c+Offset$
a :  5.7670268765345734E-01$
b :  4.9823299803773013E-01$
c :  5.4960745318930848E-01$
Offset : -2.4506729760524704E-02$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
fs(x,y);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
kill(all)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
G[0]:a0*b1*b3/b0/b2/b4$ G[1]:b0*b2*b4/b1/b3$ G[2]:b1*b3/b2/b4$ G[3]:b2*b4/b3$ G[4]:b3/b4$ G[5]:b4$
for k:4 thru 0 step -1 do print(G[k]*G[k+1])$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
Z:5 /* $ par[k]:=prod( b[2*j],j,0,k )$ imp[k]:=prod( b[2*j+1],j,0,k ) */ $
g[Z]:b[Z-1]$ b[-1]:a0$
for i:Z-1 thru 0 step -1 do print(g[i]:b[i-1]/g[i+1]);
/* [wxMaxima: input   end   ] */

/* Maxima can't load/batch files which end with a comment! */
"Created with wxMaxima"$