clf
clear
more off

% --------- BEGIN user-servicable parts ---------

% The base of the filename
BaseName='NGC7217';

% The next 3 numbers should be determined by knowledge of the true colour
% magnitudes of field stars.
Rhi=500;
Ghi=240;
Bhi=75;

Lscale=0.05;

% --------- END user-servicable parts ---------

%%% Load in the image frames
[hdr,B]=fitsload([BaseName,'-B.fits']);
[hdr,V]=fitsload([BaseName,'-V.fits']);
[hdr,R]=fitsload([BaseName,'-R.fits']);

%%% Create the RGBL image
%% The total luminosity (this should really be read from a clear band)
L=B+V+R;
% Put L on a log scale, then scale it from [0,1]
L(L<0)=L(L<0)*0+1;
L=log10(L*Lscale)/2;
L(L>1)=L(L>1)*0+1;
%% The colours
% Scale the colours to [1:256]
B_=A_scale(smooth(B,10,'gaus'),0,Bhi,256);
V_=A_scale(smooth(V,10,'gaus'),0,Ghi,256);
R_=A_scale(smooth(R,10,'gaus'),0,Rhi,256);
% The 'colour' image.  This image on its own is the standard RGB image.
Colour=double(cat(3,R_,V_,B_));
% This is the neat step.  Normalise the RGB colour for each pixel by
% the maximum value of the RGB triplet, then scale to (0,256)
% This will produce an image that has no luminosity information, just colour.
max_Colour=max(Colour,[],3);
for i=1:3, Colour(:,:,i)=256*(Colour(:,:,i)./max_Colour); end
% Multiply the colours by the luminosity
Image=cat(3,Colour(:,:,1).*L,Colour(:,:,2).*L,Colour(:,:,3).*L);
Image=uint8(Image);
Image=flipdim(Image,1);
image(Image)

imwrite(Image,[BaseName,'-LRGB.tif'],'tif')
imwrite(Image,[BaseName,'-LRGB.png'],'png')