%etamutype2.m
function [res1,res2,eta,mu]=etamutype2(m,a,b)
s=0;
M=[a b*a;b*a a];
[E,D]=eig(M);
COL=[];
for i=1:1:m
temp1=randn(1,1000000);
temp2=randn(1,1000000);
stemp=E*D.^(1/2)*[temp1;temp2];
COL=[COL stemp];
s=s+stemp(1,:).^2+stemp(2,:).^2;
end
[res1 res2]=hist(sqrt(s),100);
eta=mean(COL(1,:).*COL(2,:))/a;
h=1/(1-eta^2);
H=eta/(1-eta^2);
m=mean(s);
v=var(s);
mu=(m^2/v)*(1+(H/h)^(2));
[p1,q1,e1,m1]=etamutype2(10,2,0);
[p2,q2,e2,m2]=etamutype2(1,2,0.5);
[p3,q3,e3,m3]=etamutype2(10,1,0);
[p4,q4,e4,m4]=etamutype2(1,2,0.5);
figure
subplot(2,2,1)
plot(q1,p1)
title(strcat('eta=',num2str(e1),'mu=',num2str(m1)))
subplot(2,2,2)
plot(q2,p2,'r')
title(strcat('eta=',num2str(e2),'mu=',num2str(m2)))
subplot(2,2,3)
plot(q3,p3)
title(strcat('eta=',num2str(e3),'mu=',num2str(m3)))
subplot(2,2,4)
plot(q4,p4,'r')
title(strcat('eta=',num2str(e4),'mu=',num2str(m4)))
