clear
f=5.11;
 m=400.3;
 a=m*f/(f-1);b=m-a;
 sigma=102.66;
 lamda=m*f/((f-1)*sigma);
 d=lamda^2;
x=390:.5:410;
si=[5 47 99];
for m=1:3
for k=1:length(x),
v=si(m)^2;
w=v^2;
g=a/d;
h(k)=x(k)-b;
c(1)=2*g*w;
c(2)=-v*(v+4*h(k)*g);
c(3)=2*h(k)*(v+g*h(k));
c(4)=a^2-h(k)^2;
roots(c);
z=roots(c);
t=z(3);
n=9000;
r=h(k)-v*t;
q=a/r;
e=(-1)*h(k)*t+(1/2)*v*t^2+d*(1-q);
p=v+(r^3)/(a*d);
s(k)=(n/(2*pi*p))^1/2*exp(n*e);
end
i=min(x)
j=max(x)
l=min(s)
o=max(s)
z=trapz(x,s)
nor=1/z
s1=nor*s;
%axis([i j l o])  	
col(m,1)=plot(x,s1),
o1=nor*o
text(400,o1,'si=')


hold on
end
hold off