% eda11_03
% variance of T=1/f computed by:
%        small number approximation
%        numerical simulation
clear all;

tpi = 2*pi;

wbar = 1.3*tpi;
sigmaw = 0.2;

Nbins = 201;
wbmin = 0;
wbmax = 20;
Dwb = (wbmax-wbmin)/(Nbins-1);
wb = wbmin + Dwb*[0:Nbins-1]';

pw1 = (1/(sqrt(tpi)*sigmaw))*exp(-((wb-wbar).^2)/(2*(sigmaw^2)));

N = 10000;
wr = random('Normal',wbar,sigmaw,N,1);
pw2 = hist(wr,wb)';
pw2 = pw2 / (Dwb*sum(pw2));

clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, 10, 0, 3] );
xlabel('w (radians/s)');
ylabel('p(w)');
plot( wb, pw1, '-', 'LineWidth', 4, 'Color', [0.8,0.8,0.8] );
plot( wb, pw2, 'k:', 'LineWidth', 2 );

Tbmin = 0;
Tbmax = 2;
DTb = (Tbmax-Tbmin)/(Nbins-1);
Tb = Tbmin + DTb*[0:Nbins-1]';

Tbar = tpi/wbar;
sigmaT = sigmaw*tpi/(wbar^2);

pT1 = (1/(sqrt(tpi)*sigmaT))*exp(-((Tb-Tbar).^2)/(2*(sigmaT^2)));

Tr = tpi./wr;
pT2 = hist(Tr,Tb)';
pT2 = pT2 / (DTb*sum(pT2));

subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, 2, 0, 25] );
xlabel('T (s)');
ylabel('p(T)');
plot( Tb, pT1, '-', 'LineWidth', 4, 'Color', [0.8,0.8,0.8] );
plot( Tb, pT2, 'k:', 'LineWidth', 2 );


