% This function takes as input a pair of file names for images in "fits" format.
% The function produces an HR diagram that highlights the variable stars
% (RR-Lyrae, Cepheid, etc.).
%
% For globular cluster M15, I used this function call...
%    makeHRdiagRRlyrae( 'm15_2013-LRGB.fit', 'm15_2015-LRGB.fit', 'M15', 11, 0.5, 1.0, 'merge', 1.3); 
% For globular cluster M15, I used this function call...
%    makeHRdiagRRlyrae('m5_2005-LRGB.fit', 'm5_2006-LRGB.fit', 'M5', 11.2, 0.2, 1.0, 'merge', 1.3);

function makeHRdiagRRlyrae(fname1,fname2,objname,maxmag0,fluxthresh,colorthresh,show,z)

close all;

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

RGB2=fitsread(fname2);
R2 = RGB2(:,:,1);
G2 = RGB2(:,:,2);
B2 = RGB2(:,:,3);
L2 = (R2+G2+B2)/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');
L2convL2 = conv2(L2, Laplacian2, 'same');
R2convL2 = conv2(R2, Laplacian2, 'same');
B2convL2 = conv2(B2, Laplacian2, 'same');
thresh = 0.0006;  % 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))) ...
    & (L2convL2 > thresh * max(max(L2convL2))) ...
    & (R2convL2 > thresh * max(max(R2convL2))) & (B2convL2 > thresh * max(max(B2convL2))) ...
	;

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

figure(1);
whitebg(1,'k');
set(gcf, 'Color', [0 0 0]);
logL = log(L+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);

'hello'
[m,n] = size(k)
y = (1:m)' * ones(1,n);
x = ones(m,1) * (1:n);
size(x)
size(y)
size(k)
x = x(k);
y = y(k);

Lflux = LconvL2(k);
Rflux = RconvL2(k);
Bflux = BconvL2(k);
L2flux = L2convL2(k);
R2flux = R2convL2(k);
B2flux = B2convL2(k);

[Lfluxsort, I] = sort(Lflux);
Rfluxsort = Rflux(I);
Bfluxsort = Bflux(I);
L2fluxsort = L2flux(I);
R2fluxsort = R2flux(I);
B2fluxsort = B2flux(I);
xsort = x(I);
ysort = y(I);

min(L2fluxsort)
min(R2fluxsort)
min(B2fluxsort)

minflux = min(Lfluxsort);

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

Lfluxmag = 2.5*log10(Lfluxsort);           L2fluxmag = 2.5*log10(L2fluxsort);
maxLfluxmag = max(max(Lfluxmag));          maxL2fluxmag = max(max(L2fluxmag));
Lfluxmag = Lfluxmag-maxLfluxmag - maxmag0; L2fluxmag = L2fluxmag-maxL2fluxmag - maxmag0;
Lfluxmag = -Lfluxmag;                      L2fluxmag = -L2fluxmag;

maxmag = max(Lfluxmag);
minmag = min(Lfluxmag);
avgmag = mean(Lfluxmag(end-100:end))
avgmag2 = mean(L2fluxmag(end-100:end))
L2fluxmag = L2fluxmag - avgmag2 + avgmag;

colormag = 2.5*log10(Rfluxsort./Bfluxsort);
color2mag = 2.5*log10(R2fluxsort./B2fluxsort);
avgcol = mean(colormag(end-100:end))
avgcol2 = mean(color2mag(end-100:end))

figure(6);
plot((Lfluxmag+L2fluxmag)/2, log(abs(Lfluxmag-L2fluxmag)),'b.');
title('Randomness vs Magnitude');
xlabel('Magnitude');
ylabel('log of abs of diff of mags');
A = [ones(n,1) (Lfluxmag+L2fluxmag)/2];
b = log(abs(Lfluxmag-L2fluxmag));
ab = A \ b;
residual = sqrt(mean((b-A*ab).^2));
hold on;
plot((11:19), ab(1)+ab(2)*(11:19),'r-');
plot((11:19), ab(1)+ab(2)*(11:19) + z*residual,'r--');
drawnow;
figure(7);
plot((Lfluxmag+L2fluxmag)/2, (abs(Lfluxmag-L2fluxmag)),'b.');
title('Randomness vs Magnitude');
xlabel('Magnitude');
ylabel('abs of diff of mags');
hold on;
plot((11:0.1:19), max(exp(z*residual)*exp(ab(1)+ab(2)*(11:0.1:19)),fluxthresh),'r-');
drawnow;

hFig = figure(2);
[nn mm] = size(R);
set(hFig,'Position',[50 50 mm nn]);
whitebg(1,'k');
set(gcf, 'Color', [0 0 0]);
logR = log(R+1);              logG = log(G+1);              logB = log(B+1);
logR = logR - min(min(logR)); logG = logG - min(min(logG)); logB = logB - min(min(logB));
logR = logR+2;                logG = logG+2;                logB = logB+2;
logR = logR/max(max(logR));   logG = logG/max(max(logG));   logB = logB/max(max(logB));
Lrgb(:,:,1) = logR;           Lrgb(:,:,2) = logG;           Lrgb(:,:,3) = logB;
imagesc(Lrgb);
axis image;
axis off;
hold on;
for k = 1:length(xsort)
  color = 2.5*log10(Rfluxsort(k)/Bfluxsort(k));
  color2 = 2.5*log10(R2fluxsort(k)/B2fluxsort(k)) - avgcol2 + avgcol;
  if abs(Lfluxmag(k) - L2fluxmag(k)) > fluxthresh ...
     && log(abs(Lfluxmag(k) - L2fluxmag(k))) ...
	- (ab(1) + ab(2)*(Lfluxmag(k)+L2fluxmag(k))/2) > z*residual ...
     && abs(color - color2) < colorthresh
	plot(xsort(k),ysort(k),'ro');
  end
end
hold off;
drawnow;

hFig = figure(3);
[nn mm] = size(R);
set(hFig,'Position',[50 50 mm nn]);
whitebg(1,'k');
set(gcf, 'Color', [0 0 0]);
logR2 = log(R2+1);               logG2 = log(G2+1);               logB2 = log(B2+1);
logR2 = logR2 - min(min(logR2)); logG2 = logG2 - min(min(logG2)); logB2 = logB2 - min(min(logB2));
logR2 = logR2+2;                 logG2 = logG2+2;                 logB2 = logB2+2;
logR2 = logR2/max(max(logR2));   logG2 = logG2/max(max(logG2));   logB2 = logB2/max(max(logB2));
Lrgb2(:,:,1) = logR2;            Lrgb2(:,:,2) = logG2;            Lrgb2(:,:,3) = logB2;
imagesc(Lrgb2);
axis image;
axis off;
hold on;
for k = 1:length(xsort)
  color = 2.5*log10(Rfluxsort(k)/Bfluxsort(k));
  color2 = 2.5*log10(R2fluxsort(k)/B2fluxsort(k)) - avgcol2 + avgcol;
  if abs(Lfluxmag(k) - L2fluxmag(k)) > fluxthresh ...
     && log(abs(Lfluxmag(k) - L2fluxmag(k))) ...
	- (ab(1) + ab(2)*(Lfluxmag(k)+L2fluxmag(k))/2) > z*residual ...
     && abs(color - color2) < colorthresh
	plot(xsort(k),ysort(k),'ro');
  end
end
hold off;
drawnow;

hFig = figure(5);
set(hFig,'Position',[10 10 960 831]);
%set(hFig,'Position',[10 10 768 665]);
axis ij;
whitebg(5,'k');
set(gcf, 'Color', [0 0 0]);
hold on;
for k = 1:n
        bright = sqrt((Lfluxmag(k)-maxmag)/(minmag-maxmag));

        maxmag = max(maxmag, 18);
        diam = ((Lfluxmag(k)-maxmag)/(12-maxmag));
	color = 2.5*log10(Rfluxsort(k)/Bfluxsort(k));
	color2 = 2.5*log10(R2fluxsort(k)/B2fluxsort(k)) - avgcol2 + avgcol;
        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 strcmp(show , 'merge') == 1
	  if abs(Lfluxmag(k) - L2fluxmag(k)) > fluxthresh ...
	     && log(abs(Lfluxmag(k) - L2fluxmag(k))) ...
	        - (ab(1) + ab(2)*(Lfluxmag(k)+L2fluxmag(k))/2) > z*residual ...
	     && abs(color - color2) < colorthresh
          plot([color color2], [Lfluxmag(k) L2fluxmag(k)], 'o-' ...
		    ,'Color',[r g b] ...
		    ,'MarkerFaceColor',[r g b] ...
		    ,'MarkerSize',marksize ...
		    ,'MarkerEdgeColor',[r g b] ...
    		    );
	  else
          plot((color+color2)/2, (Lfluxmag(k)+L2fluxmag(k))/2, 'o-' ...
		    ,'MarkerFaceColor',[r g b] ...
		    ,'MarkerSize',marksize ...
		    ,'MarkerEdgeColor',[r g b] ...
    		    );
	  end
	elseif strcmp(show , 'nolines') == 1
	  if abs(Lfluxmag(k) - L2fluxmag(k)) > fluxthresh ...
	     && log(abs(Lfluxmag(k) - L2fluxmag(k))) ...
	        - (ab(1) + ab(2)*(Lfluxmag(k)+L2fluxmag(k))/2) > z*residual ...
	     && abs(color - color2) < colorthresh
          plot([color color2], [Lfluxmag(k) L2fluxmag(k)], 'o' ...
		    ,'MarkerFaceColor',[r g b] ...
		    ,'MarkerSize',marksize ...
		    ,'MarkerEdgeColor',[r g b] ...
    		    );
	  else
          plot((color+color2)/2, (Lfluxmag(k)+L2fluxmag(k))/2, 'o-' ...
		    ,'MarkerFaceColor',[r g b] ...
		    ,'MarkerSize',marksize ...
		    ,'MarkerEdgeColor',[r g b] ...
    		    );
	  end
	else
          plot([color], [Lfluxmag(k)], 'o-' ...
		    ,'MarkerFaceColor',[r g b] ...
		    ,'MarkerSize',marksize ...
		    ,'MarkerEdgeColor',[r g b] ...
    		    );
          plot([color2], [L2fluxmag(k)], 'o-' ...
		    ,'MarkerFaceColor',[b r g/2] ...
		    ,'MarkerSize',marksize ...
		    ,'MarkerEdgeColor',[b r g/2] ...
    		    );
	end
end
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{white} Color index (B-R)');
ylabel('\color{white} Brightness (apparent magnitude)');
title(['\color{white} Hertzsprung-Russell Diagram for ' objname]);
hold off;

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