function photo = imread_fits(FITSfile, frame)
% photo = imread_fits(FITSfile, frame)
% load a single frame from an astronomical FITS file,
% such as the ones at
%   http://www.telescope.org/
% .
% The first frame is "frame 0".
% The actual frame loaded is modulo how many frames there are;
% so calling
%   photo = imread_fits('H:\download\telescope\ex_00801.fit', -1);
%   views(photo);
% actually loads the *last* frame.
% FITSfile must be the *full* path+filename of a ".fits" file.
% FIXME: this should really be merged into the IMREAD
% routine, to automatically detect FITS files.
% WARNING: still experimental.
% WARNING: this trims the rightmost column off the images,
% to get rid of the (instrument failure ?)
% bogus zeros there
% on some of the http://www.telescope.org/ files
% (in particular, these files:
%  ex_01201.fit ring nebula
%  ex_00801.fit galaxy
% have this problem)
% See also IMREAD, VIEWS, IMREAD_PCT.

% change log:
% 1999-01-21:DAV: started.
% David Cary <d.cary at ieee.org>

% open the file
% if the images look funny, omit the ieee-be flag

fp = fopen(FITSfile,'r','ieee-be');

s = char(fread(fp, 80, 'char')')

if(fp==-1)
    error(['Can''t open the file ' FTSfile ' in imread_fit.m']);
else
    if(~strncmp('SIMPLE  = ', s, 10))
        error([...
                '??? Error using ==> imread_fit\n'...
                'Unable to determine the file format\n'...
                '['...
                FITSfile...
                '] Not in FITS format.'...
            ])
    else
        naxis3 = 1; % default of 1 image
        done = 0;
        % FITS file header
        % is always a multiple of 36 records of 80 bytes each.
        while(~done)
            for i = 1:36,
                s = char(fread(fp, 80, 'char')');
                % disp(s)
                
                if(strncmp('BITPIX  = ', s, 10))
                    [bits_per_pixel] = sscanf(s(10:80), '%i')
                end
                
                if(strncmp('NAXIS   = ', s, 10))
                    [naxis] = sscanf(s(10:80), '%i')
                end
                
                if(strncmp('NAXIS1  = ', s, 10))
                    [naxis1] = sscanf(s(10:80), '%i')
                end
                
                if(strncmp('NAXIS2  = ', s, 10))
                    [naxis2] = sscanf(s(10:80), '%i')
                end
                
                if(strncmp('NAXIS3  = ', s, 10))
                    [naxis3] = sscanf(s(10:80), '%i')
                end
                
                if(or(...
                        strncmp('END             ', s, 10),...
                        feof(fp)...
                        )),
                    done = 1;
                end
            end
        end
        
        
        % s = char(fread(fp, 80, 'char')');
        % disp(['[' s ']'])
        % s = char(fread(fp, 80, 'char')');
        % disp(['[' s ']'])
        % datastart = ftell(fp);
        
        if( ~isequal( 16,bits_per_pixel ) )
            bits_per_pixel
            error(...
               'Sorry -- reading this type of image not yet implemented. Check (d.cary at ieee) for a later version.');
         end
        if( ~isequal( 1,naxis3 ) )
            naxis3
            disp('multiples image file')
        end
        
        photo=zeros(naxis1, naxis2);
        
        % assumes 16 bits per pixel
        % offset = datastart + (naxis1*naxis2*2).*mod(frame, naxis3);
        % fseek(fp,offset,'bof');
        
        % No transpose needed to make it look right on screen.
        % fseek(fp, datastart, 'bof');
        photo = fread(fp,[naxis1 (naxis2-1)],'uint16');
        % !!! Throws away the last column of data !!!
        
        
    end;
    fclose(fp);
    
end

% end imread_fits.m