function makeHRdiagram(fname,objname,maxmag0)

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');
RconvL2 = conv2(R, Laplacian2, 'same');
BconvL2 = conv2(B, 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))) ...
    & (RconvL2 > thresh * max(max(RconvL2))) & (BconvL2 > thresh * max(max(BconvL2))) ...
	;

% 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);
Rflux = RconvL2(k);
Bflux = BconvL2(k);

[fluxsort, I] = sort(flux);
Rfluxsort = Rflux(I);
Bfluxsort = Bflux(I);

minflux = min(fluxsort);

[n,m] = size(flux);
color = 2.5*log10(Rfluxsort./Bfluxsort);

fluxmag = 2.5*log10(fluxsort);
maxfluxmag = max(max(fluxmag));
fluxmag = fluxmag-maxfluxmag - maxmag0;
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;
cumflux = 0;
cumRflux = 0;
cumBflux = 0;
m = 0;
printflag = 1;
for k = 1:n
    cumflux = cumflux + fluxsort(k);
    cumRflux = cumRflux + Rfluxsort(k);
    cumBflux = cumBflux + Bfluxsort(k);
    m = m+1;
    if (cumflux > random('Poisson',10)*minflux) 
	avgflux = cumflux/m;
	avgRflux = cumRflux/m;
	avgBflux = cumBflux/m;
	avgfluxmag = 2.5*log10(avgflux);
	avgfluxmag = avgfluxmag - maxfluxmag - maxmag0;
	avgfluxmag = -avgfluxmag;
        bright = sqrt((avgfluxmag-maxmag)/(minmag-maxmag));
%	bright = 1;
        maxmag = max(maxmag, 18);
        diam = ((avgfluxmag-maxmag)/(12-maxmag));
	color = 2.5*log10(avgRflux/avgBflux);
        r = max(min(1 + 0.7*color,1),0);
        b = max(min(1 - 0.7*color,1),0);
        g = 1 - 0.5*(1-min(r,b));
        r = bright*r;
        g = bright*g;
        b = bright*b;
        marksize = (0.0001+diam*2);
	if m == 1 && printflag == 1
	    objname
	    avgfluxmag
	    printflag = 0;
	end
        plot(color, avgfluxmag, 'o' ...
		    ,'MarkerFaceColor',[r g b] ...
		    ,'MarkerSize',marksize ...
		    ,'MarkerEdgeColor',[r g b] ...
    		    );
	cumflux = 0;
	cumRflux = 0;
	cumBflux = 0;
	m = 0;
    end
end
%if strcmp(objname , 'M29') == 1
%axis([-1.0 2.0 8 17]);         % for M29
%elseif strcmp(objname , 'M92') == 1
%axis([-0.8 1.0 10 20]);         % for M92
%elseif strcmp(objname , 'M13') == 1
%axis([-1.2 1.4 11 19]);         % for M13
%else
if strcmp(objname , 'M35 and NGC2158') == 1
axis([-2 2 8 17]);         % for M35
elseif strcmp(objname , 'NGC7789') == 1
axis([-2 2 8 17]);         % for M35
else
axis([-2 2 10 21]);         % for M35
end
xlabel('Color index (B-R)');
ylabel('Brightness (apparent magnitude)');
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);
