clear
clf

% ------- START of user input --------

% List of filters used.  Choose from C U B V R & I.
Filters='BVR';

% The exposures corresponding to the set of flats for each filter.
Exposure=[120 10 1.6];

% The location of the flats
FlatDir='../Flats/';

% The number of flats for each filter
Nflats=[3 3 3];

% The location of the darks
DarkDir='../Flats/';

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



%-------- END of user input ( generally ) -----------
NColours=length(Nflats);
for colour=1:NColours

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

 % Read the files in and subtract the dark
 for i=1:Nflats(colour);
  File=[FlatDir,'Flat-',Filters(colour),'-',num2str(Exposure(colour)),'s-',int2str(i),'.fits'];
  [hdr,imgs(:,:,i)]=fitsload(File);
  imgs(:,:,i)=(imgs(:,:,i)-dark);
 end

 % Balance the images for variations in exposure due to things
 % like variation in the voltage to the lamps
 [N,M]=size(imgs(:,:,1));
 for i=1:Nflats(colour)
  [n,x]=hist(reshape(imgs(:,:,i),N*M,1),1000);
  n_max_pos=find(n==max(n));
  span_n=[n_max_pos-10:n_max_pos+10];
  P=gauss_fit(x(span_n),n(span_n),[integrate(n)*(x(2)-x(1)),std(x(span_n)),x(n_max_pos)]);
  imgs(:,:,i)=imgs(:,:,i)/P(3);
 end
 % The previous can sometimes be done just as well using:
 % for i=1:Nflats(colour)
 % imgs(:,:,i)=imgs(:,:,i)/mean(mean(imgs(:,:,i)));
 % end
 

 % Get the median
 img=median(imgs,3);

 imagesc(img,[0.98 1.035]);
 axis equal

 % Output the flat to FITS
 name=['Flat_',Filters(colour),'.fits'];
 hdr=fitsremovekey(hdr,'BITPIX');
 hdr=fitsaddkey(hdr,'BITPIX',-32,2,'4 byte float');
 C=clock;
 hdr=fitsaddkey(hdr,'COMMENT',[],0,[int2str(Nflats(colour)),...
                ' flats combined to create this flat.']);
 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 colour