clear
pack

t100=clock;

eta=.667;h=6.625*10^(-34);c=3*10^8;l=1.173*10^(-6);
ni=c/l;p=1.0*10^(-9.1);
g=80;kef=.04;t1=(22/9)*10^(-9);
q=1.6*10^(-19);
fe=kef*g+(2-1/g)*(1-kef);
ne=p*eta*t1/(h*ni);
sigma=(g*q*sqrt(ne*fe))/t1;me=(ne*g*q)/t1;
lambda=(sqrt(ne*fe))/(fe-1);
a1=[4.5 21 50 80 82 60 30 9 3.18];
a=(1/339.68).*a1;
%m1=me.*a;
%sigma1=sigma.*a;
%m1(1)=me;
%sigma1(1)=sigma;
x=(0:4:800);
mi=me*10^9
sigmai=sigma*10^9
0ui=1+(x-mi)./(lambda*sigmai);
wi=(1./(sqrt(2*pi)*sigmai.*ui.^1.5)).*exp(-1.*(x-mi).^2 ./(2*sigmai^2 .*ui));
y=x(1:2:200);
cc=wi(2:2:200);
pn=8.*cc;tt=sum(pn)
r=cc(51:100);
v2=ones(1,50);
d=cc(50).*v2;v1=cc(1:50);
pm=8.*[d r];
v11=(1/cc(50)).*v1;
v3=[v11 v2];ss=sum(pm)
s1=cumsum(pm);
t11=clock
nl=25;nc=9;
ru=ss.*rand(nl,nc);
t1a=clock
et=zeros(nl,nc);es=zeros(nl,nc);pon=zeros(nl,nc);sz=[0];
fr=[sz s1];n1=length(fr)-1;contor=0;di=zeros(1,nl);
for k=1:n1

for m=1:nl
%di=all(et');
if ((di(m)>=1)&(m<nl))
m=m+1;
elseif(m<=nl)
for n=1:nc
if ((et(m,n)>=1)&(n<nc))
n=n+1;
else
aa=fr(k);bb=fr(k+1);
if(ru(m,n)>=aa)&(ru(m,n)<bb)
et(m,n)=1;
es(m,n)=y(k)+(ru(m,n)-fr(k))/pm(k);
pon(m,n)=v3(k);
else
contor=contor+1;
end
end
end
end
end
di=all(et');
end

b1=all (et)
durata=etime(clock,t11) 
ys=(22/9)*es*a';
wgt=cumprod(pon');




h1=[64 0 0 0 0 0 0 0];
h2=[64 32 0 0 0 0 0 0];
h3=[0 32 0 0 0 0 0 0];
