disp('Problem 9')

function y=I(alpha); y=p(alpha).*q(alpha); end
function dy=dI(alpha); dy=p(alpha).*dq(alpha)+dp(alpha).*q(alpha); end

function y=p(alpha); y=2+sin(10*alpha); end
function dy=dp(alpha); dy=10*cos(10*alpha); end

function y=q(alpha);
global method alfa
alfa=alpha;
switch method
  case 1, y=mprule('f1',0,2);
  case 2, y=quad('f2',0,2);
  case 3, y=fourint('f3',alpha/2,alpha/(2*pi));
  case 4, y=p9G(alpha);
end
end

function dy=dq(alpha);
global method alfa
alfa=alpha;
switch method
  case 1, dy=mprule('df1',0,2);
  case 2, dy=quad('df2',0,2);
  case 3, dy1=fourint('df3',alpha/2,alpha/(2*pi));
          dy2=fourint('df4',alpha/2,alpha/(2*pi)+0.5);
          dy=dy1+dy2;
  case 4, dy=dG(alpha);
end
end

% int(a:b) func(x) fx
function y=mprule(func,a,b)
global neval
h=(b-a)/neval; y=h*sum(feval(func,a+h*(0.5:neval)));
end

% int(0:inf) func(x) sin(omega*x + pi*theta) dx
function [y,ya]=fourint(func,omega,theta)
global bigM
[t,w]=q_ossinf(bigM,omega,theta); 
f=feval(func,t);
y=sum(w.*f);
ya=sum(w.*abs(f))*eps;
end

function y=f1(x)
global alfa
y=x.^alfa.*sin(alfa./(2-x));
end

function y=df1(x)
global alfa
y=log(x).*f1(x)+x.^alfa.*cos(alfa./(2-x))./(2-x);
end

function y=f2(t)
global alfa
z=t+i*t.*(2-t); dz=1+2i*(1-t);
y=real(-i*exp(alfa*(log(z)+i./(2-z))).*dz);
end

function y=f3(t)
global alfa
x=2*t./(1+t); dx=2./(1+t).^2;
y=x.^alfa.*dx;
end

function y=df2(t)
global alfa
z=t+i*t.*(2-t); dz=1+2i*(1-t);
a=(log(z)+i./(2-z));
y=real(-i*a.*exp(alfa*a).*dz);
end

function y=df3(t)
global alfa
x=2*t./(1+t); dx=2./(1+t).^2;
y=log(x).*x.^alfa.*dx;
end

function y=df4(t)
global alfa
x=2*t./(1+t); dx=2./(1+t).^2;
y=x.^alfa./(1+t);
end

function [alf,f]=solve(mth)
global method; method=mth; set_options;
alf(1)=pi/4-0.005; alf(2)=pi/4;
f(1)=dI(alf(1)); f(2)=dI(alf(2));
for k=3:20,
  alf(k)=alf(k-1)-f(k-1)*(alf(k-1)-alf(k-2))/(f(k-1)-f(k-2));
  f(k)=dI(alf(k));
  if abs(f(k))>=abs(f(k-1)), break; end
end
end

% symmetry axis of parabola through three given points is at x(3)+h
function [h,th]=extreme(x,y)
h=diff(x); d=diff(y)./h; th=0.5*d(2)/(d(1)-d(2));
h=th*h(1)+(th-0.5)*h(2);
end

% Maximize I(alpha), using specified method.
% Do maxk iterations (default: stop at best value) 
function [alf,f]=maximize(mth,maxk)
global method; method=mth; set_options;
if nargin<2, maxk=20; end
alf=pi/4-[0.01; 0.005; 0];  f=alf; for j=1:3, f(j)=I(alf(j)); end
for k=4:maxk,
  j=k-[3 2 1]; h=extreme(alf(j),f(j)); 
  alfnew=alf(k-1)+h; fnew=I(alfnew);
  if nargin<2, if fnew <= f(k-1), break, end; end 
  alf(k)=alfnew; f(k)=fnew;
end
end

function set_options
global bigM 
bigM=16; quad_options('abs',1e-10); quad_options('rel',0);
end

function [alf,f]=findzero(mth)
global method; method=mth; set_options;
alf=pi/4-[0.005; 0]; f=alf; for j=1:2, f(j)=dI(alf(j)); end
for k=2:20,
  h=-(alf(k)-alf(k-1))*f(k)/(f(k)-f(k-1)); 
  alfnew=alf(k)+h; fnew=dI(alfnew);
  if abs(fnew) >= abs(f(k)), break, end;
  alf(k+1)=alfnew; f(k+1)=fnew;  
end
end

function y=dG(alpha,h); if nargin<2, h=0.001; end
dh=(p9G(alpha+h)-p9G(alpha-h))/(2*h);
d2h=(p9G(alpha+2*h)-p9G(alpha-2*h))/(4*h);
y=dh+(dh-d2h)/3;
end