function makeHRdiagram(fname,objname,maxmag)

close all;

RGB=fitsread(fname);
R = RGB(:,:,1);
G = RGB(:,:,2);
B = RGB(:,:,3);
L = (R+G+B)/3;

m1  = [  0 1 -1 ];
m2  = [ -1 1  0 ];
m3  = [  0 1 -1 ]';
m4  = [ -1 1  0 ]';
m5  = [ -1 0 0 ;
         0 1 0 ;
	 0 0 0 ]';
m6  = [ 0 0  0 ;
        0 1  0 ;
	0 0 -1 ]';
m7  = [  0 0 0 ;
         0 1 0 ;
	-1 0 0 ]';
m8  = [ 0 0 -1 ;
        0 1  0 ;
	0 0  0 ]';
m9  = [ 0 0 0 1 -1 ];
m10  = [ -1 1 0 0 0 ];
m11  = [  0 0 0 1 -1 ]';
m12  = [ -1 1 0 0  0 ]';
Laplacian  = [ 0 -1 0;  -1 4 -1; 0 -1 0 ];
Laplacian2 = [ 0 -1 -1 -1  0 ;  
              -1  0  2  0 -1 ; 
              -1  2  4  2 -1 ; 
              -1  0  2  0 -1 ; 
               0 -1 -1 -1  0 ]; 
Lconv1 = conv2(L, m1, 'same');
Lconv2 = conv2(L, m2, 'same');
Lconv3 = conv2(L, m3, 'same');
Lconv4 = conv2(L, m4, 'same');
Lconv5 = conv2(L, m5, 'same');
Lconv6 = conv2(L, m6, 'same');
Lconv7 = conv2(L, m7, 'same');
Lconv8 = conv2(L, m8, 'same');
Lconv9 = conv2(L, m9, 'same');
Lconv10 = conv2(L, m10, 'same');
Lconv11 = conv2(L, m11, 'same');
Lconv12 = conv2(L, m12, 'same');
LconvL = conv2(L, Laplacian, 'same');
LconvL2 = conv2(L, Laplacian2, 'same');
thresh = 0.0003;  % 0.001 for M29
k = (Lconv1 > 0) ...
  & (Lconv2 > 0) ...
  & (Lconv3 > 0) ...
  & (Lconv4 > 0) ...
  & (Lconv5 > 0) ...
  & (Lconv6 > 0) ...
  & (Lconv7 > 0) ...
  & (Lconv8 > 0) ...
  & (Lconv9 > 0) ...
  & (Lconv10 > 0) ...
  & (Lconv11 > 0) ...
  & (Lconv12 > 0) ...
    & (LconvL > thresh * max(max(LconvL))) & (LconvL2 > thresh * max(max(LconvL2))) ...
	;

% trim off 2 pixels around the edge
k(:,1:2) = 0;
k(:,end-1:end) = 0;
k(1:2,:) = 0;
k(end-1:end,:) = 0;

figure(1);
whitebg(1,'k');
set(gcf, 'Color', [0 0 0]);
logL = log(L+1);
logL = logL - min(min(logL));
%logL = 65535*logL/max(max(logL));
%logL = log(logL+1);
%logL = logL - min(min(logL));
logL = logL/max(max(logL));
Lrgb(:,:,1) = logL - logL.*k;
Lrgb(:,:,2) = logL + (1-logL).*k;
Lrgb(:,:,3) = logL - logL.*k;
imagesc(Lrgb);
axis image;
axis off;
set(gcf, 'InvertHardCopy', 'off');
fname=['HRdiag/tickmarks' objname];
print('-f1','-dpng', '-r50', fname);
print('-f1','-depsc', fname);

flux = LconvL2(k);
color = log(R(k)./B(k));
%figure(4);
%axis xy;
%whitebg(4,'w');
%semilogy(color, flux, 'k.');
%%axis([-1.1 1.6 100 500000]);
%xlabel('Color index (B-R)');
%ylabel('Relative flux (adu)');
%title('Hertzsprung-Russell Diagram for '+objname);

fluxmag = 2.5*log10(flux);
fluxmag = fluxmag-max(max(fluxmag)) - maxmag;
fluxmag = -fluxmag;
maxmag = max(fluxmag);
minmag = min(fluxmag);
figure(5);
axis ij;
%h = gcf;
%set(h,'Color',[0 0 0]);
whitebg(5,'k');
set(gcf, 'Color', [0 0 0]);
hold on;
[n, m] = size(color);
for k = 1:n
    diam = ((fluxmag(k)-maxmag)/(minmag-maxmag));
    bright = (diam);
    maxmag = max(maxmag, 18);
    diam = ((fluxmag(k)-maxmag)/(12-maxmag));
    r = max(min(1 + 0.7*color(k),1),0);
    b = max(min(1 - 0.7*color(k),1),0);
    g = 1 - 0.5*(1-min(r,b));
    r = bright*r;
    g = bright*g;
    b = bright*b;
    marksize = (0.0001+diam*2);
    plot(color(k), fluxmag(k), 'o' ...
		    ,'MarkerFaceColor',[r g b] ...
		    ,'MarkerSize',marksize ...
		    ,'MarkerEdgeColor',[r g b] ...
    		    );
end
%axis([-1.0 2.0 8 17]);         % for M29
%axis([-0.8 1.0 10 20]);         % for M92
xlabel('Color index (B-R)');
ylabel('Brightness (mag)');
title(['\color{white} Hertzsprung-Russell Diagram for ' objname]);
hold off;

set(gcf, 'InvertHardCopy', 'off');
fname=['HRdiag/HRdiag' objname];
print('-f5','-dpng', '-r100', fname);
print('-f5','-depsc', fname);
%print('-f5','-dpdf', fname);