global A, global B, global epsr
B=1.36e-03;
epsr=1.0e-06;

kk=[];zz=[];A0=[];
for A=[0.21:0.001:0.25 0.26:0.0025:1.96 1.965:0.001:2]
 k0=0;zz0=[];
 if (A<0.25)
 [z0 tz0 k0 zz0]=rk4st("fsurf",[0;0;0;0;0],0,0.015,80000);
 else
 [z0 tz0 k0 zz0]=rk4st("fsurf",[0;0;0;0;0],0,0.003/A,80000);
 endif
 kk=[kk k0];
 A0=[A0 A];
 zz=[zz zz0];
endfor

hh=loglog(zz(4,:).^(1/3),zz(2,:),zz(9,:).^(1/3),zz(7,:),zz(14,:).^(1/3),zz(12,:),zz(19,:).^(1/3),zz(17,:),zz(24,:).^(1/3),zz(22,:),zz(29,:).^(1/3),zz(27,:),zz(34,:).^(1/3),zz(32,:),zz(39,:).^(1/3),zz(37,:),zz(44,:).^(1/3),zz(42,:),(pi*zz(2,:)/6.*(3+zz(2,:))).^(1/3),zz(2,:),":b",(pi*zz(7,:)/6.*(3+zz(7,:))).^(1/3),zz(7,:),":b")
set(hh,"linewidth",2)
axis([1e-1 30 3e-1 1e2])
daux=0.01:0.01:100;

hh=loglog(zz(4,:).^(1/3),zz(2,:),zz(9,:).^(1/3),zz(7,:),zz(14,:).^(1/3),zz(12,:),zz(19,:).^(1/3),zz(17,:),zz(24,:).^(1/3),zz(22,:),zz(29,:).^(1/3),zz(27,:),zz(34,:).^(1/3),zz(32,:),zz(39,:).^(1/3),zz(37,:),zz(44,:).^(1/3),zz(42,:),(pi*daux/6.*(3+daux)).^(1/3),daux,":k");
set(hh,"linewidth",2)
axis([5e-1 30 3e-1 1e2])