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