clear
clf

%-------- Start of user input ---------------

% For referencing the data
BaseName='NGC7217';
BaseDir='../NGC7217/';

% The filters used.  Choice of C, U, B, V, R, & I
Filters='BVR';

% Where the darks are
DarkDir='../Darks/';

% The exposures of the darks and image frames for each colour
Exposure=[200 200 200];

% The number of darks for each colour
Ndarks=[3 3 3];

% The number of image frames
Nfiles=[3 3 3];

% The position from which to extract an image
% Remember that x is the second index!
StarFramex=[130:160];
StarFramey=[295:325];

%-------- End of user input ( generally ) -----------

NColours=length(Filters);

for colour=1:NColours

 %%%Read in the darks and median
 for i=1:Ndarks
 File=[DarkDir,'Dark-',num2str(Exposure(colour)),'s-',int2str(i),'.fits'];
  [hdr,darks(:,:,i)]=fitsload(File);
 end
 dark=median(darks,3);

 % Read in the Flat
 [hdr_flat,flat]=fitsload(['Flat-',Filters(colour),'.fits']);

 % Read the files in, subtract the dark, and flat-field
 for i=1:Nfiles;
  filename=[BaseDir,BaseName,'-',num2str(Exposure(colour)),'s-',...
            Filters(colour),'-',int2str(i),'.fits'];
  [hdr,imgs(:,:,i)]=fitsload(filename);
  imgs(:,:,i)=(imgs(:,:,i)-dark)./flat;
 end

 %%% Register the images with each other
 % Centre of frame
 meanx=mean(StarFramex);
 meany=mean(StarFramey);
 for i=1:Nfiles
  % Find x and y by fitting a 2-D Gaussian to the selected star.
  % Warning messages are often returned.
  [bias,sigma,ampl,x,y]=fit_gaussian(imgs(StarFramey,StarFramex,i)*1000);
  % Find the shift to apply
  dx=meanx-(x+StarFramex(1)-1);
  dy=meany-(y+StarFramey(1)-1);
  % Apply the shift
  imgs(:,:,i)=shift(imgs(:,:,i),dy,1);
  imgs(:,:,i)=shift(imgs(:,:,i),dx,2);
 end

 %%% Subtracts the approximate sky BG from all the images.
 %% Fit the sky to each colour with a 2nd-order polynomial plane.
 %% Then subtract the 'sky' from each colour.
 sigmas=2;
 sigma=zeros(Nfiles(colour),1);
 for i=1:Nfiles
  B=imgs(:,:,i);
  mask=MaskStars(B,sigmas);
  [Coef,Fit]=FitSurface(B,mask,2);
  imgs(:,:,i)=B-Fit;
  % Store the sigma for the sky for later use.
  sigma(i)=std(B(find(mask)));
 end

 %%% Balance the images for variations in sky transparency
 %%% Find the typical scale and normalise by that.
 %%% To find the typical scale, use the first image as a baseline, and then
 %%% find mean(img(:,:,i)./img(:,:,1)).  To eliminate the noise from dominating,
 %%% use only those parts of the data that are > 10 sigma for each
 % Initialise some arrays
 [N,M]=size(imgs(:,:,1));
 Ps=zeros(Nfiles(colour),1);
 % Main loop
 img_base=squeeze(imgs(:,:,1));
 Ps(1)=1;
 for i=2:Nfiles
  % Extract one image at a time
  tmp_img=squeeze(imgs(:,:,i));
  % Find those pixels > 10 sigma above the sky
  signal=find(tmp_img>10*sigma(i)); % find (cts-<BG>) > 10*sigma(BG)
  % Note, if a cosmic-ray pixel is selected by the signal criterion, then it
  % will have a very small value in the base, so the ratio will come out as 0.
  ratios=img_base(signal)./tmp_img(signal);
  % clip the data before finding the mean to elliminate the small values.
  scale=mean(clip(ratios,3)); 
  imgs(:,:,i)=scale*tmp_img;
  Ps(i)=scale;
 end

 %%% Get the median
 %% Most of the previous work was done so this step works!
 img=median(imgs,3);

 %%% Output the image to FITS
 name=[BaseName,'-',Filters(colour),'.fits'];
 hdr=fitsremovekey(hdr,'BITPIX');
 hdr=fitsaddkey(hdr,'BITPIX',-32,2,'4 byte float');
 C=clock;
 hdr=fitsaddkey(hdr,'COMMENT',[],0,[int2str(Nfiles),...
                ' images combined to create this.']);
 hdr=fitsaddkey(hdr,'COMMENT',[],0,['Created ',int2str(C(1)),...
                '-',int2str(C(2)),'-',int2str(C(3))]);
 hdr=fitsremovekey(hdr,'BZERO');
 hdr=fitsremovekey(hdr,'BSCALE');
 hdr=fitsremovekey(hdr,'PLATESCA');
 hdr=fitsremovekey(hdr,'RA');
 hdr=fitsremovekey(hdr,'DEC');
 hdr=fitsremovekey(hdr,'EQUINOX');
 if(exist(name))
  eval(['!rm -f ',name])
 end
 status = fitswrite(hdr,img,name);

end % loop over all colours