Plotting echelle spectra

[admin]

Hello,


echelle spectra, when merged, are long high resolution spectra. In order to fully use the screen or paper when potting them, I have writtent this small code to cut the spectrum into N sections.

It will be at the end part of my 'otz_plot' generic function (-e option) but until then, here it is as you may find this useful...

Please do not hesitate to reply here if you see ways to improve this code, I'm a true beginner in MatLab!


Cordialement,
Olivier

Code: Select all

% otz_plot_echelle.m
%
% Display and create a PDF of a graph of an echelle spectrum
%
% v1.0 / 20160816 - published in ARAS forum
%
% 

clear all % efface toute les variables

% Parameters
NbGraphs = 10;  % number of graphs (ie: the spectrum is cut in 'NbGraphs' parts)
MarginPercent = 0; % [0-50] spectrum margin to add in each subgraphs (in % of the individual graph size)

% Here is the FITS file which includes the spectrum & key information
% One option is to ask the user to point toward the FITS file
[FileName, PathName] = uigetfile ('*.fits'); 
FileWay = [PathName, FileName];
% Antoher option is to write here the FITS file name including directory
%FileWay = 'C:\Users\Olivier\Documents\MATLAB\Demo_FITS\_vvcep_20160719_956_full.fits';

% Extract FITS header and key infos from there
FITSHeader = fitsinfo(FileWay);
Data = fitsread(FileWay);
% Look for the size of the header (nb of keywords)
EnteteMax = max(size(FITSHeader.PrimaryData.Keywords));
% Extract some data from the header
% (to be improved, see rfitsinfo.m work in progress)
for l=1:EnteteMax
    k = FITSHeader.PrimaryData.Keywords{l,1};
    
    switch k
        case 'NAXIS1'
            NbPoints = int32(FITSHeader.PrimaryData.Keywords {l,2});
        case 'CRVAL1'
            Origine = FITSHeader.PrimaryData.Keywords {l,2};
        case 'CDELT1'
            Pas = FITSHeader.PrimaryData.Keywords {l,2};
        case 'OBJNAME'
            Objet = FITSHeader.PrimaryData.Keywords {l,2};
        case 'DATE-OBS'
            Date = FITSHeader.PrimaryData.Keywords {l,2};
        case 'OBSERVER'
            Observateur = FITSHeader.PrimaryData.Keywords {l,2};
        case 'BSS-VHEL'
            VHel = FITSHeader.PrimaryData.Keywords {l,2};
        case 'EXPTIME'
            ExpTimeTotal = FITSHeader.PrimaryData.Keywords {l,2};
        case 'EXPTIME2'
            ExpTime = FITSHeader.PrimaryData.Keywords {l,2};
        case 'BSS_SITE'
            SiteObs = FITSHeader.PrimaryData.Keywords {l,2};
        case 'BSS_INST'
            Instrument = FITSHeader.PrimaryData.Keywords {l,2};
        case 'JD-OBS'
            JDobs = FITSHeader.PrimaryData.Keywords {l,2};
        case 'JD-MID'
            JDmid = FITSHeader.PrimaryData.Keywords {l,2};
    end
end


% Define the individual graph size (in Nb of points) and the margin
NbPointsPerGraph = idivide(NbPoints, int32(NbGraphs), 'round');
if MarginPercent <= 0
    NbPointsMargin = 0;
elseif MarginPercent > 50
    NbPointsMargin = idivide(NbPointsPerGraph, int32(2), 'round');
else    
    NbPointsMargin = idivide(NbPointsPerGraph, int32(100/MarginPercent), 'round');
end

% Open a full screen figure
figure('units','normalized','outerposition',[0 0 1 1])

% For each graph...
for k = 1 : NbGraphs
    
    aMin = max(1,((k-1)*NbPointsPerGraph)-NbPointsMargin);
    aMax = min(NbPoints,(k*NbPointsPerGraph)+NbPointsMargin);

    clear x;
    clear y;
    
    % Création des données x et y du spectres
    for a = aMin : aMax
        Lambda = Origine + double(a-1) * Pas;
        x (a-aMin+1) = Lambda;
        y (a-aMin+1) = Data (1,a);
    end
    LambdaMin = Origine + (aMin-1) * Pas;
    LambdaMax = Origine + (aMax-1) * Pas;
    
    % Affiche le graphe
    ax=subplot(NbGraphs,1,k);
    plot(x, y, 'r');

    % Set no margin on the X axis display
    ax.XLim = [-inf inf];
    ax.XAxis.FontSize = 6;
    ax.XAxis.TickLength = [0.0020 0.0050];
    ax.YLim = [-inf inf];
    ax.YAxis.FontSize = 6;
    ax.YAxis.TickLength = [0.0020 0.0050];
    
    % Display X axis label below last graph
    xlhand = get(gca,'XLabel');
    if k == NbGraphs
        set(xlhand,'string','Wavelength (A)','fontsize', 12);
    end
    
    % Display title above the first graph
    if k == 1
        Titre2 = [Objet, ' (', Date, ' / ', ExpTime, ') @ ', SiteObs, ' w/ ', Instrument,'; (c) ', Observateur];
        title(Titre2, 'Interpreter','none');
        %xlabel('Wavelenght in Â');
        %ylabel('Relative flux');
    end
end

% sauve le graphe en pdf
orient landscape;
print ('graph', '-painters', '-dpdf', '-bestfit');
print ('graph', '-painters', '-depsc');
print ('graph', '-painters', '-dsvg');