function f = rhoplot(rfile,zfile,rhofile,ipnt,rmax,outfile)
%----------------------------------------
% Plot la densité de nucléons dans le plan (Oxz)
% rfile et zfile = nom des fichiers contenant les coordonnées r (en fait x) et z
% rhofile = non du fichier contenant les densités
% ipnt = choix de la densité à ploter (1=neutron, 2=proton, 3=tot, 4=diff)
% rmax = bornes max de la densité interpolée
% outfile = nom du fichier de sortie
%-----------------------------------------
% Lecture des data (coordonnees)
r = textread(rfile,'%n');
z = textread(zfile,'%n');
% Definition des bornes dans les directions z et r
ngh = size(z,1);
ngl = size(r,1);
% Lecture des data : a et b non utilisés, rhoin densité totale de test, puis densité neutron, proton,totale
[a,b,rhon,rhop,rhot,rhodif] = textread(rhofile,'%n%n%n%n%n%n');
% On met les variables sous forme de matrices
rhon = reshape(rhon,[ngh,ngl])';
rhop = reshape(rhop,[ngh,ngl])';
rhot = reshape(rhot,[ngh,ngl])';
rhodif = reshape(rhodif,[ngh,ngl])';
% Intervales sur lesquels seront interpolées les densités
rr = [ -rmax : 0.025 : rmax ];
zz = [ -rmax : 0.025 : rmax ];
%Definition du mesh d'interpolation
[zzz,rrr] = meshgrid( zz,rr );
% Interpolation des densités sur le mesh
rhon = interp2( z, r, rhon, zzz, rrr, 'spline' );
rhop = interp2( z, r, rhop, zzz, rrr, 'spline' );
rhot = interp2( z, r, rhot, zzz, rrr, 'spline' );
rhodif = interp2( z, r, rhodif, zzz, rrr, 'spline' );
% Choix de la quantité à ploter
if(ipnt==1)
rhoplot = rhon
end
if(ipnt==2)
rhoplot = rhop
end
if(ipnt==3)
rhoplot = rhot
end
if(ipnt==4)
rhoplot = rhodif
end
% Declaration de la figure et des axes
f = figure('position',[100 0 1000 1000],'paperpositionmode','auto', 'color','w','InvertHardcopy','off');
PLotHandle = axes('Position',[.24 .1 .75 .8],'Visible','off','FontSize',25);
% On trace la densité
pcolor(rrr,zzz,rhoplot)
%Generation de la colormap
shading flat
colorDepth = 200;
colormap(jet(colorDepth))
%On impose aux axes d'etre égaux en longueur
axis square
%Définition des labels
xlabel('x(fm)')
ylabel('z(fm)')
% On sauvegarde
hold on
%On ajoute des contours d'isodensité
ContourHandle = contour(rrr,zzz,rhoplot,'k--');
%On met une échelle de couleur
ch = colorbar('location','East','FontSize',15); 
set(ch, 'YColor',[1 1 0.6])
% On incruste un plot 3D pour le style
Plot3DHandle = axes('Position',[.01 .60 .40 .40],'Visible','off'); %le nom des axes
surf(rrr,zzz,rhoplot)  %la fig 3D
shading flat % option de rendu
set(Plot3DHandle,'Box','on')   %on affiche la boite 3D 
% Pour des raisons de lisibilité, on ne veut ni le nom des axes, ni les valeurs correpondantes
set(Plot3DHandle,'XTick',[])
set(Plot3DHandle,'XTickLabel',[])
set(Plot3DHandle,'YTick',[])
set(Plot3DHandle,'YTickLabel',[])
set(Plot3DHandle,'ZTick',[])
set(Plot3DHandle,'ZTickLabel',[])
% On rentre à la main les bornes de x et y
set(Plot3DHandle,'XLim',[-rmax rmax])
set(Plot3DHandle,'YLim',[-rmax rmax])
%set(PeaksPlotHandle,'ZLim',[0 0.22])
% On imprime ca dans un fichier

%print(f,outfile,'-depsc','-painters')  % format vectoriel mais ne gere pas la transparence 
%print(f,outfile,'-depsc','-zbuffer','-r200')
print(f,outfile,'-djpeg','-zbuffer','-r200')
%print(f,'-depsc2','-painters',outfile)
