mu=0.50; psi=0.33;
v0=0:0.001:1;
v1=transpose(v0);
h(:,1)=mu*v1(:,1).*(1-v1(:,1)).^(0.95);
h(:,2)=-1*h(:,1);
m0(:,1)=psi*v1(:,1).*(1-v1(:,1)).^(0.95);
m0(:,2)=-1*m0(:,1);
vv0(:,1)=v1;
vv0(:,2)=v1;
hmax=0;
mmax=0;
for i=1:1:length(h)
if h(i,1)>hmax
hmax=h(i,1);
poshmax=(i);
end
if m0(i,1)>mmax
mmax=m0(i,1);
posmmax=(i);
end
end
tol=0.0015*hmax;
h1=1/4*hmax;
j=0;
for i=1:1:length(h)
if abs(h(i,1)-h1)<tol
j=j+1;
hh1=h(i,1);
poshh1(j,1)=(i);
end
end
h2=1/2*hmax;
j=0;
for i=1:1:length(h)
if abs(h(i,1)-h2)<tol
j=j+1;
hh2=h(i,1);
poshh2(j,1)=(i);
end
end
h3=4/5*hmax;
j=0;
for i=1:1:length(h)
if abs(h(i,1)-h3)<tol
j=j+1;
hh3=h(i,1);
poshh3(j,1)=(i);
end
end
h4=1/3*hmax;
j=0;
for i=1:1:length(h)
if abs(h(i,1)-h4)<tol
j=j+1;
hh4=h(i,1);
poshh4(j,1)=(i);
end
end
h5=2/3*hmax;
j=0;
for i=1:1:length(h)
if abs(h(i,1)-h5)<tol
j=j+1;
hh5=h(i,1);
poshh5(j,1)=(i);
end
end
vh1=zeros(101,2);
vh2=zeros(101,2);
vh3=zeros(101,2);
vh4=zeros(101,2);
vh5=zeros(101,2);
vh1(:,1)=h1;
vh1(:,2)=h1;
vh2(:,1)=h2;
vh2(:,2)=h2;
vh3(:,1)=h3;
vh3(:,2)=h3;
vh4(:,1)=h4;
vh4(:,2)=h4;
vh5(:,1)=h5;
vh5(:,2)=h5;
vv1=transpose(linspace(v1(poshh1(1,1)),v1(poshh1(2,1)),101));
vv2=transpose(linspace(v1(poshh2(1,1)),v1(poshh2(2,1)),101));
vv3=transpose(linspace(v1(poshh3(1,1)),v1(poshh3(2,1)),101));
vv4=transpose(linspace(v1(poshh4(1,1)),v1(poshh4(2,1)),101));
vv5=transpose(linspace(v1(poshh5(1,1)),v1(poshh5(2,1)),101));
h1v=transpose(linspace(h1,h1,101));
h2v=transpose(linspace(h2,h2,101));
h3v=transpose(linspace(h3,h3,101));
h4v=transpose(linspace(h4,h4,101));
h5v=transpose(linspace(h5,h5,101));
m1(:,1)=sqrt(psi^2.*vv1(:,1).^2.*(1-vv1(:,1)).^(2*0.95)-psi^2.*(h1v(:,1)./mu).^2);
m2(:,1)=sqrt(psi^2.*vv2(:,1).^2.*(1-vv2(:,1)).^(2*0.95)-psi^2.*(h2v(:,1)./mu).^2);
m3(:,1)=sqrt(psi^2.*vv3(:,1).^2.*(1-vv3(:,1)).^(2*0.95)-psi^2.*(h3v(:,1)./mu).^2);
m4(:,1)=sqrt(psi^2.*vv4(:,1).^2.*(1-vv4(:,1)).^(2*0.95)-psi^2.*(h4v(:,1)./mu).^2);
m5(:,1)=sqrt(psi^2.*vv5(:,1).^2.*(1-vv5(:,1)).^(2*0.95)-psi^2.*(h5v(:,1)./mu).^2);
m1(:,2)=-1.*m1(:,1);
m2(:,2)=-1.*m2(:,1);
m3(:,2)=-1.*m3(:,1);
m4(:,2)=-1.*m4(:,1);
m5(:,2)=-1.*m5(:,1);
m1(1,:)=0;
m1(101,:)=0;
m2(1,:)=0;
m2(101,:)=0;
m3(1,:)=0;
m3(101,:)=0;
m4(1,:)=0;
m4(101,:)=0;
m5(1,:)=0;
m5(101,:)=0;
vv1(:,2)=vv1(:,1);
vv2(:,2)=vv2(:,1);
vv3(:,2)=vv3(:,1);
vv4(:,2)=vv4(:,1);
vv5(:,2)=vv5(:,1);
figure(2);
plot(real(m1),vv1,'Linestyle',':','Linewidth',2,'Color',[0.4660 0.6740 0.1880])
hold on
plot(real(m4),vv4,'Linestyle','--','Linewidth',1,'Color',[0.4660 0.6740 0.1880])
hold on
plot(real(m2),vv2,'Linestyle','-.','Linewidth',1,'Color',[0.4660 0.6740 0.1880])
hold on
plot(real(m5),vv5,'Linestyle','--','Linewidth',2,'Color',[0.4660 0.6740 0.1880])
hold on
plot(real(m3),vv3,'Linestyle',':','Linewidth',1,'Color',[0.4660 0.6740 0.1880])
hold on
plot(m0,vv0,'Linestyle','-','Linewidth',2,'Color','b')
axis([-0.2 0.2 0.0 1.0])
xlabel('M/BV_{max}');
y=ylabel('V/V_{max}','Rotation',0);
set(y, 'Units', 'Normalized', 'Position', [-0.12, 0.5, 0]);
legend('H=1/4Hmax','H=1/3Hmax','H=1/2Hmax','H=2/3Hmax','H=4/5Hmax','H=0','Location','northeast');
legend('boxoff');
