function f = rhoplotJPawesome(rfile,zfile,rhofile,outfile)

r = textread(rfile,'%n');
z = textread(zfile,'%n');

ngh = size(z,1);
ngl = size(r,1);

[a,b,rhoin] = textread(rhofile,'%n%n%n');

rho = reshape(rhoin,[ngh,ngl])';

%if( allquadrants ~= 0 )
%   r  = cat( 1, - flipud(r), r );
%   z  = cat( 1, - flipud(z), z );
%   r1 = cat( 2, fliplr(rho), rho );
%   r2 = cat( 1, flipud(r1), r1 );
%   rhonew = r2;
%end
%rr = [ -rmax : 0.025 : rmax ];
%zz = [ -rmax : 0.025 : rmax ];
%nr = 2*rmax/0.025 + 1
nr = ngl
[zzz,rrr] = meshgrid( z,r );

%rhonew = interp2( z, r, rhonew, zzz, rrr, 'spline' );
rhojp = rho((nr+1)/2,:)
rhojp2 = rho(:,(nr+1)/2)
%rhojp = rhonew(:,(nr+1)/2)
%f = figure();
f = figure('position',[100 0 1000 1000],'paperpositionmode','auto', 'color','w','InvertHardcopy','off');
%PeaksCPlotHandle = axes('Visible','off','FontSize',30);
PeaksCPlotHandle = axes('Position',[.24 .1 .75 .8],'Visible','off','FontSize',25);

% Make the contour plot
pcolor(rrr,zzz,rho)
%pcolor(peaks)
shading flat
colorDepth = 200;
colormap(jet(colorDepth))

set(gca, 'XTick', [-10 -8 -6 -4 -2 0 2 4 6 8 10]);
set(gca, 'YTick', [-10 -8 -6 -4 -2 0 2 4 6 8 10]);

axis square


xlabel('x(fm)')
ylabel('z(fm)')

hold on
ContourHandle = contour(rrr,zzz,rho,'k--');
%clabel(ContourHandle,'FontSize',15,'Color','w','Rotation',0)
%text_handle = clabel(ContourHandle,'FontSize',15,'Color',[.796 0.663 .529],'Rotation',0);
%set(text_handle,'BackgroundColor',[1 1 .6],'Edgecolor',[.7 .7 .7])
ch = colorbar('location','East','FontSize',15); 
%caxis([0.0 0.19])
%set(ch, 'ytick', (0:0.01:0.19),'YColor',[1 1 0.6])
set(ch, 'YColor',[1 1 1.0])
cutx = 0;
cuty = 0;
%peakdata = [rrr,zzz,rhonew];
%plot(cutx*ones(1,nr),r,'Color',[0.545 0.373 0.373],'LineWidth',2)

%plot(zz,cuty*ones(1,nr),'k')



%plot(cutx*ones(1,length(zzz)),[1:length(zzz)],'y')
%plot([1:length(rrr)],cuty*ones(1,length(rrr)),'r')
%RightAxesHandle = axes('Position',[.820 .155 .12 .8]);
%RightAxesHandle = axes('Position',[.855 .1 .12 .8],'XTick',[0:0.1:0.2]);
%[RightAxesHandle] = plot(rhojp,zz,'Color','m');
%plot(rhojp,z,'Color',[0.545 0.373 0.373],'LineWidth',2);
%[RightAxesHandle, CutXHandle, CutYHandle] = plotyy(peakdata(:,cutx+1),[1:length(peakdata)], peakdata(cuty+1,:),[1:length(peakdata)]);
%set(RightAxesHandle(1),'XLim',[0 0.22],'FontSize',15);
%xlabel('\rho(fm^{-3})')


%set(CutXHandle,'Color','y')
%set(CutYHandle,'Color','r')
%set(RightAxesHandle(1),'YTick',[],'YTickLabel',[],'YColor','y')
%set(RightAxesHandle(2),'YTick',[10 30 40],'YTickLabel', ...
%['a'; 'b'; 'c'],'YColor','r')
%set(RightAxesHandle(1),'XTick',[],'XTickLabel',[])
%set(RightAxesHandle(2),'XTick',[],'XTickLabel',[])
PeaksPlotHandle = axes('Position',[.01 .60 .40 .40],'Visible','off');

surf(rrr,zzz,rho)
shading flat
%ch = colorbar('location','East','FontSize',15); 
%caxis([0 0.19])
%set(ch, 'ytick', (0:0.02:0.12),'YColor',[1 1 0.6])
set(PeaksPlotHandle,'Box','on')
set(PeaksPlotHandle,'XTick',[])
set(PeaksPlotHandle,'XTickLabel',[])
set(PeaksPlotHandle,'YTick',[])
set(PeaksPlotHandle,'YTickLabel',[])
set(PeaksPlotHandle,'ZTick',[])
set(PeaksPlotHandle,'ZTickLabel',[])
set(PeaksPlotHandle,'XLim',[-10 10])
set(PeaksPlotHandle,'YLim',[-10 10])
set(PeaksPlotHandle,'ZLim',[0 1.00])
%set(RightAxesHandle, 'ytick', (0:0.1:0.2))
%print -dpsc2 outfile.eps
%print(f,outfile,'-depsc','-painters')
%print(f,outfile,'-depsc','-zbuffer','-r200')
print(f,outfile,'-djpeg','-zbuffer','-r200')
%print(f,'-depsc2','-painters',outfile)
