批量读取文件夹下的DR8光谱数据，并完成显示保存图片以及统计Z和Obj的分布情况

function [Z Obj] = read_dr8_dir(path,option,show_n, save_img, hist_info)
%Author:shizhixin
%Email:szhixin@gmail.com
%Blog:http://blog.csdn.net/shizhixin
%函数功能：批量读取文件夹下的DR8光谱数据，
%并完成显示保存图片以及统计Z和Obj的分布情况


% option='dered'
% path = 'dr8\'
% show_n=-1;
% save_img = 1;
% hist_info = 1;
% [Z Obj] = read_dr8_dir(path,option,show_n,save_img,hist_info);
%
%参数说明
% path: 需要读取的文件夹PATH
% option: 是否退红移动，'dered' or 'undered'
% show_n: 处理文件夹文件的个数，方便处理文件夹中部分文件；
% if show_n=-1 then 处理所有文件
% if show_n>总数，then 处理所有文件
% save_img:
% 是否显示并保存数据的图于当前路径下，取1 or 0
% hist_info:
% 是否获取头文件中的Z Obj,并做直方图，取1 or 0；
% if hist_info == 1,[Z Obj]返回其对应的值，否则[Z Obj]都为0
% Obj的值表示含义：
% GALAXY  Obj = 1;
% QSO Obj = 2;
% STAR Obj=3;
% else Obj=4;
                
Z = zeros(show_n,1);%红移
Obj = zeros(show_n,1);%类型


files = dir([path '*.fits']);
n = length(files)
if show_n > n
    disp 'input show_n is larger than the number of files!set show_n = n!'
    show_n = n;
end
if show_n == -1
    show_n=n;
end


for i = 1:show_n
    i
    filename = [path files(i).name]
    [hed lamd flux z] = read_dr8_flux(filename,option);
    z
    %-----------------plot line and save image----------------
    if save_img == 1
        plot(lamd,flux);


        hold on;
        plot_dot_line(flux);


        fits_name = files(i).name;
        save_name = strtok(files(i).name,'.');
        print(gcf,'-dpng',[path save_name '.png']);
        hold off
        title(save_name);
    end% if save_img
    %---------------------------------------
    if hist_info == 1
        [Z(i),Obj(i)] = get_hed_info(hed);
    end%if hist_info
    %-------------------------------
end%for


if hist_info == 1
    hist_fit_info(Obj, Z);
end
end%fun




%------------------------------------------------------------------------
function plot_dot_line(flux)
min_flux = min(flux);
max_flux = max(flux);
line_set = [4862 4959 5007 6564 6583 6716 6731];%[H-BETA,OIII, OIII, H-ALPHA, NII, SII,SII]
num_line = size(line_set,2);
for plot_i=1:num_line
    plot([line_set(plot_i):line_set(plot_i)],[min_flux:5:max_flux],'r-');
    hold on
end
end
%------------------------------------------------------------------------


%------------------------------------------------------------------------
function [z,obj] = get_hed_info(hed)
size_hed = size(hed,1);
for j=1:size_hed
    switch hed(j,1:8)
        
        case 'RAOBJ   '
            RAOBJ=strtrim(hed(j,11:30));


        case 'DECOBJ  '
            DECOBJ=strtrim(hed(j,11:30));


        case 'QUALITY '
            QUALITY=strtrim(hed(j,11:30));


        case 'OBJTYPE '
            OBJTYPE=strtrim(hed(j,12:19));
            if strcmp(OBJTYPE, 'GALAXY')
                obj = 1;
            elseif strcmp(OBJTYPE, 'QSO')
                obj = 2;
            elseif strcmp(OBJTYPE, 'STAR')
                obj=3;
            else
                obj=4;
            end


        case 'Z       '
            red=strtrim(hed(j,11:30));
            z = str2double(red(2:end-1));
    end%switch
end%for
end
%------------------------------------------------------------------------
%------------------------------------------------------------------------
function hist_fit_info(Obj, Z)
figure;
hist(Obj, [1:4]);
title('The hist of ObjType');
figure;
Z(find(Z>5))=0;
hist(Z, [0:0.05:1]);
title('The hist of Z');
end
%------------------------------------------------------------------------

function [hed,lamd,flux,z]=read_dr8_flux(filename,option)
%Author:shizhixin
%Email:szhixin@gmail.com
%Blog:http://blog.csdn.net/shizhixin


% filename = 'dr8/spec-2951-54592-0433.fits';
% option = 'undered'  or 'dered'


% DR8 FITS文件
% 文件头1   2880字节 END结束
% 文件头2   N*2880字节 END结束
% 数据区：
% flux，double，长度为NAXIS2
% best_fit double 长度为NAXIS2
% wavelength double 长度为NAXIS2
% and_mask int32 长度为NAXIS2
% or_mask int32 长度为NAXIS2
% inverse_variance double长度为NAXIS2
% 不是按块，而是按顺序依次存储. 
% flux1, best_fit1....flux2, best_fit2...




fid=fopen(filename,'r','ieee-be');
fseek(fid, 2880, -1); %跳过文件开头2880个字节


%计算文件头中包含的2880的块数
read = true;
while read == true
    fseek(fid, 80, 0);
    string = fread(fid,3,'*char');
    fseek(fid, -3, 0);
    pos = ftell(fid);
    if string(1)=='E' & string(2)=='N' & string(3)=='D'
        read = false;
    end
end
pos = pos+3;%'D'的后面一个位置
chunk_num = ceil(pos/2880);%计算文件头中包含的2880的块数，包括第一个2880


fseek(fid, 2880, -1); %跳过文件开头2880个字节
hed_unit_num = (2880*(chunk_num-1))/80; %读取文件头，文件头总大小为5*2880，每个单元80个字节
hed=fread(fid,[80,hed_unit_num],'*char')';%读取出来的数据存入数组, in column order


for i=1: hed_unit_num
    switch hed(i,1:8)
        case 'NAXIS2  '
            dim=str2num(hed(i,10:30));%dim;number of pixel
        case 'COEFF0  '
            start=str2num(hed(i,12:19));%start,wavelentgh of first pixe


        case 'COEFF1  '
            deat=str2num(hed(i,12:19));%deat,dispersion per pixel


        case 'CRVAL1  '
            start=str2num(hed(i,12:19));
        case 'CD1_1   '
            deat = str2num(hed(i,12:19));
        case 'Z       '
            z=str2num(hed(i,12:19));%redshift
    end
end


lamd0=start:deat:start+deat*(dim-1);
lamd0 = 10.^lamd0';


%-----read flux of spectrum-----
Oset = 2880*chunk_num;
status=fseek(fid,Oset,-1);
for i = 1:dim
    flux(i) = fread(fid,1,'double');
    best_fit(i) = fread(fid,1,'double');
    lamd(i) = fread(fid,1,'double'); % 和lamd0一样
    and_mask(i) = fread(fid,1,'int32');
    or_mask(i) = fread(fid,1,'int32');
    inverse_variance(i) = fread(fid,1,'double');
end
fclose(fid);


switch option
    case 'undered'
        lamd=lamd;
    case 'dered'
        %---dered--------
        lamd=lamd./(1+z);
end


%  plot(lamd, flux);
%  figure
%  plot(lamd, best_fit);


end %function