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^(-8.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;m=(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=m.*a;
%sigma1=sigma.*a;
%m1(1)=m;
%sigma1(1)=sigma;
x=linspace(0,800,256);
mi=m*10^9
sigmai=sigma*10^9
ui=1+(x-mi)./(lambda*sigmai);
wi=(1./(sqrt(2*pi)*sigmai.*ui.^1.5)).*exp(-1.*(x-mi).^2 ./(2*sigmai^2 .*ui));
y=[0 x(1:2:256)];
cc=wi(2:2:256);
pn=8.*cc;tt=sum(pn)
r=cc(65:128);
v2=ones(1,64);
d=wi(119).*v2;v1=cc(1:64);
pm=8.*[d r];
v11=(1/wi(119)).*v1;
v3=[v11 v2];ss=sum(pm)
s1=cumsum(pm);
t11=clock
nl=100;nc=9;fr=[s1];
ru=ss.*rand(nl,nc);
t1a=clock
es=zeros(nl,nc);pon=zeros(nl,nc);
for j=1:nl 
  for h=1:nc  nrz=ru(j,h);
      zm=1;zM=128;zme=64;
         for i=1:7
            if nrz<=s1(zme) zM=zme;zme=fix((zm+zM)/2);
            else zm=zme;zme=fix((zm+zM)/2);
            end
         end
         es(j,h)=y(zm)+(ru(j,h)-fr(zm))/pm(zm); 
         pon(j,h)=v3(zm);                                                    
  end
end
durata=etime(clock,t11) 
pack
ys=(22/9)*es*a';
wgt=prod(pon');
pondh=zeros(1,128);
ysM=max(ys);ysm=min(ys);lint=(ysM-ysm)/127;
for l=1:length(ys)
 zm=0;zM=128;zme=64;nrl=ys(l);
         for i=1:7
            if nrl<=zme*lint zM=zme;zme=(zm+zM)/2;
            else zm=zme;zme=(zm+zM)/2;
            end
         end
         pondh(zM)=pondh(zM)+wgt(l);
end
probe=pondh./nl;
ril=1:length(probe);
plot(ril,probe)  
 