% to calculate phase integral int(pdq) for
% harmonic oscillator for hydrogen
global cay mu energy
	cay=0.37;mu=918.55;
clear action
clear e
for ie=1:50
	energy=5.e-4+ie*.001;
   % initialize initial guess for turning point (only for first energy,
   % otherwise the initial guess is the turning point at the previous
   % energy
   if ie==1
      rturn=0.5
   end
	rturn=fzero('vquadratic',rturn); 
	s(ie)=2*quad('pquadratic',-rturn+2e-8,rturn-2e-8);
    if ie==1
		disp([cay mu energy ie rturn s(ie)])
    end
	e(ie)=energy;
end
s=s/(2*pi);

[ex mn]=min(abs(s-0.5));
cc=polyfit(e(mn-1:mn+1),s(mn-1:mn+1)-.5,2);
ee(1)=min(roots(cc));

[ex mn]=min(abs(s-1.5));
cc=polyfit(e(mn-1:mn+1),s(mn-1:mn+1)-1.5,2);
ee(2)=min(roots(cc));


