function G=Gnew(v,vf,vs,vm)
% Gnew=vg(v) for phase transition. v is the volume fraction and g(v) is the
% radial distribution function.
% usage:: G=Gnew(v,[vf],[vs],[vm])
% 
% vf is the volume fraction at the fluidus (default vf=0.65). 
% vs is the volume fraction at the solidus (default vs=0.72). 
% vm is the maximum volume fraction (default vm=pi/sqrt(3)/2). 

% revision history:
% 12/29/09 Mark D. Shattuck <mds> Gnew.m 

if(~exist('vf','var'))
  vf=.65;
end
if(~exist('vs','var'))
  vs=.72;
end
if(~exist('vm','var'))
  vm=pi/sqrt(3)/2;
end

Gf=vf*(16 - 7*vf)./( 16*(1 - vf).^2);  % G at vf
Sf=(8+vf)/8/(1-vf)^3;                  % G' at vf
Ss=Gf/(vm-vs);                         % G' at vs

a=(Sf+Ss)/(vf-vs)^2;                   % make above true at transision
vb=(Ss*vf+Sf*vs)/(Sf+Ss);               

G=a*(v-vf).*(v-vb).*(v-vs)+Gf;         % G in transition

ii=v<vf;
G(ii)=v(ii).*(16-7*v(ii))./(1-v(ii)).^2/16;   % G below

ii=v>vs;
G(ii)=Gf*(vs-vm)./(v(ii)-vm);          % G above


