timour_matlab/mscripts/frmax.m

function frm=frmax(p,le,r)
if (p~=1.)
a=p-1.;
b=(p+1.-2.*p.*le).*r./le;
c=(p.*le+1.-2.*p).*r./le;
d=(p.*le-1.).*r.*r./le./le;
h=(3.*a.*c-b.*b)./(9.*a.*a);
k=(2.*b.*b.*b-9.*a.*b.*c+27.*a.*a.*d)./(27.*a.*a.*a);
th=(acos(-k./(2.*(-h).^1.5)))./3.;
end %if

if (p>1)
  frm=2.*sqrt(-h).*cos(th+4*pi/3)-b./(3.*a);
elseif (p<1)
  frm=2.*sqrt(-h).*cos(th)-b./(3.*a);
elseif (p==1)
  frm=0.5.*sqrt(0.25+2.*r./le)+0.25;
elseif (p==0)
  frm=sqrt(r./le)
end
return
1	function frm=frmax(p,le,r)
2	if (p~=1.)
3	a=p-1.;
4	b=(p+1.-2.p.le).*r./le;
5	c=(p.le+1.-2.p).*r./le;
6	d=(p.le-1.).r.*r./le./le;
7	h=(3.a.c-b.b)./(9.a.*a);
8	k=(2.b.b.b-9.a.b.c+27.a.a.d)./(27.a.a.a);
9	th=(acos(-k./(2.*(-h).^1.5)))./3.;
10	end %if
11
12	if (p>1)
13	frm=2.sqrt(-h).cos(th+4pi/3)-b./(3.a);
14	elseif (p<1)
15	frm=2.sqrt(-h).cos(th)-b./(3.*a);
16	elseif (p==1)
17	frm=0.5.sqrt(0.25+2.r./le)+0.25;
18	elseif (p==0)
19	frm=sqrt(r./le)
20	end
21	return