%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Instable Fermi-Bose mixture [Chui04], numerical integration.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;						addpath C:\Users\Philippe\Physik\Appendix\Appendix;		Constant;
figure(1);					clf;
col='krgbmcy';
m =    7*uA;																				%mass
as=-27.3*aB;				g=4*pi*hP^2*as/m;												%boson-boson scattering length
wR=2*pi*50;					a =sqrt(hP/m/wR);												%trap
N =[500:500:3000];	
r =[0:.01:1]*50e-6;																			%atom numbers
w =[.05:.05:10]*a;																			%variad radius of condensate
for jj=1:length(N)
	for j=1:length(w)
		U   =m/2*wR^2*r.^2;																	%harmonic potential
		phi = sqrt(N(jj)*1/w(j)^3/pi^1.5)*exp(-r.^2/2/w(j)^2);								%wave function

		Heff=hP^2/2/m* (-r/w(j)^2.*phi).^2 + U.*phi.^2 + g/2*phi.^4;	
		EB(j)=sum(sum(4*pi*r.^2.* Heff)) *(r(2)-r(1));

% 		subplot(221),plot(r*1e6,U/hP/wR)
% 		set(gca,'FontSize',14);
% 		xlabel('\fontsize{14}{\itr}  (\mum)');
% 		ylabel('\fontsize{14}{\itU/h\omega}');
% 		shading flat;
% 		subplot(222),plot(r*1e6,phi.^2*1e-6)
% 		set(gca,'FontSize',14);
% 		xlabel('\fontsize{14}{\itr}  (\mum)');
% 		ylabel('\fontsize{14}{\itn_B}   (cm^{-3})');
% 		shading flat;
% 		pause(.01);
	end
	semilogx(w/a,EB/hP/wR,col(jj),'LineWidth',1.4)
	set(gca,'xLim',[.1 10],'yLim',[0 40000],'FontSize',16,'LineWidth',1.4);
	xlabel('\fontsize{18}{\it\sigma_r/a_{trap}}');
	ylabel('\fontsize{18}{\itE} [{\it\psi}_0] / {\ith\nu}');
	hold on;
end

