function [k,a,yksikohta] = kertoimet(alpha,R)
%Funktio laskee kertoimet
%
%a_k=alpha(1-alpha)gamma(k+alpha+1)/k!gamma(alpha+1)
%
%arvoilla k=0,1,...,R vektoriin a ja plottaa niistä kuvan.
%Funktio laskee myös arvion luvulle k, jolle a_k=1.
%
%Kutsu funktiota komennolla
%
%[k,a,yksikohta] = kertoimet(alpha,R);
%
%esimerkiksi
%
%[k,a,yksikohta] = kertoimet(0.5,10^3);

k=0:R;
a=alpha*(1-alpha)*(k.*beta(k,alpha+1)).^(-1);
b=a.^(1/alpha); %Funktion b(k) kuvaaja lähestyy suoraa L.
figure
hold on
plot(k,a);
plot(k,b,'r');

title('$$a_k=\alpha(1-\alpha)\frac{\Gamma(k+\alpha+1)}{k!\Gamma(\alpha+1)}$$','interpreter','latex','fontsize',20)
xlabel('$$k$$','interpreter','latex','fontsize',20)
ylabel('$$a_k$$','interpreter','latex','fontsize',20)
h=legend('$a_k$','$a_k^{\frac{1}{\alpha}}$','Location','NorthWest');
set(h,'Interpreter','latex')
set(h,'FontSize',20)

r=round(R/2);
kk=(b(R)-b(r))/(R-r);
%
%Suoran L yhtälö on y=kk*(x-R)+b(R), joten kohta, jossa L leikkaa suoran
%y=1, voidaan nyt ratkaista ja saadaan arvio luvulle k, jolla a_k=1.
%Huomaa, että b(k)=1 jos ja vain jos a(k)=1.

yksikohta=(1-b(R))/kk+R-1;
Yk=ceil(yksikohta);
yk=floor(yksikohta);

aYk=alpha*(1-alpha)*(Yk.*beta(Yk,alpha+1)).^(-1);
ayk=alpha*(1-alpha)*(yk.*beta(yk,alpha+1)).^(-1);

disp(['yksikohta=' num2str(yksikohta,10)])
disp(['a(' num2str(Yk) ')=' num2str(aYk,10)])
disp(['a(' num2str(yk) ')=' num2str(ayk,10)])
end

