function [output_image, bzero, bscale] = fitsread(filename,n_hdu)
%FITSREAD reads a FITS image into a Matlab variable. 
%
%  [data, bzero, bscale]=fitsread(filename)
%  [data, bzero, bscale]=fitsread(filename,h_hdu)
%
% The variable 'n_hdu' explicity tells FITSREAD how many
% sets of 36 header 'cards' precede the data. One normally
% lets the function figure this out on its own.
%
% Final data values are to be computed by data*bscale+bzero
%
% Known problems
%	Will only deal with 2D fits files.
%
% Version 2.0 
% R. Abraham, Institute of Astronomy, Cambridge University
%
% Modified R.G.Lane 18 November 1997 to use strcmp 
%                               to return bzero, bscale
% EEE University of Canterbury                              

%The first few cards of the first set of 36 cards must give information 
%in a pre-defined order (eg. the data format, number of axes, 
%size etc) but after that the header keywords can come in any
%order. The end of the last card giving information is flagged by 
%the keyword END, after which blank lines are added to pad the 
%last set of cards so it also contains 36 cards. After the last card
%the data begins, in a format specified by the BITPIX keyword.
%The dimensions of the data are specified by the NAXIS1 and
%NAXIS2 keywords.
%
%Reference:   NASA/Science Office of Standards and Technology
%           "Definition of the Flexible Image Transport System"
%                            NOST 100-1.0
%
%           This and other FITS documents are available on-line at:
%           <A HREF="http://www.gsfc.nasa.gov/astro/fits/basics_info.html">http://www.gsfc.nasa.gov/astro/fits/basics_info.html</A>

%Set flag indicating unknown number of HDUs
if nargin<2
	n_hdu='unknown';
end;

if isstr(n_hdu)
	n_hdu=upper(n_hdu);
end;

%Allow user to specify a few keywords to indicate number of card sets
if strcmp(n_hdu,'HST') | strcmp(n_hdu,'STSCI') | strcmp(n_hdu,'MDS')
   n_hdu=6;
end
if strcmp(n_hdu,'FREI')
   n_hdu=2;
end

%Open the file
fid=-1;
if ~isstr(filename)
	filename=setstr(filename);
end;
if (isempty(findstr(filename,'.'))==1)
	filename=[filename,'.fits'];
end
[file,message] = fopen(filename,'r','l');
if file == -1
	error(message);
end


%First five cards must contain specific information
[d,simple,d]=parse_card(setstr(fread(file,80,'uchar')'));
[d,bitpix,d]=parse_card(setstr(fread(file,80,'uchar')'));
[d,naxis,d]=parse_card(setstr(fread(file,80,'uchar')'));
[d,naxis1,d]=parse_card(setstr(fread(file,80,'uchar')'));
[d,naxis2,d]=parse_card(setstr(fread(file,80,'uchar')'));


if n_hdu=='UNKNOWN'
	%Keep reading cards until one turns up with the keyword 'END'.
	n_card=5;
	keyword=' ';
	while(~strcmp(upper(deblank(keyword)),'END'))
      		n_card=n_card+1;
    		card=setstr(fread(file,80,'uchar')');
                keyword=parse_card(card);
        	if (strcmp(upper(deblank(keyword)),'BZERO'))
			[keyword,bzero,comment]=parse_card(card);
			end
        	if (strcmp(upper(deblank(keyword)),'BSCALE'))
    		        [keyword,bscale,comment]=parse_card(card);
			end
	end;
	%Go past the blank lines of padding before the start of the data
	n_blanks = 36 - rem(n_card,36);
	dummy=fread(file,n_blanks*80,'uchar');
else
	dummy=fread(file,((n_hdu*36)-5)*80,'uchar'); 
end;

%Read the data. Note big-endian switch is used (Mac and UNIX!
if bitpix==-64
   X=fread(file,naxis1*naxis2,'float32','b');
elseif bitpix==-32
   X=fread(file,naxis1*naxis2,'float','b');
elseif bitpix==8
   X=fread(file,naxis1*naxis2,'uint8','b');   
elseif bitpix==16
   X=fread(file,naxis1*naxis2,'short','b');   
elseif bitpix==32
   X=fread(file,naxis1*naxis2,'long','b');      
else
   error('data type specified by BITPIX keyword is not -64, -32, 8, 16, or 32');
end;  

%Clean up and output data 
fclose(file);
output_image=rot90(reshape(X,naxis1,naxis2));