Plotting spectral profile

Discussion et échanges sur l'utilisation de MatLab en spectroscopie
Post Reply
Olivier Thizy
Posts: 370
Joined: Sat Sep 24, 2011 10:52 am
Location: in the french Alps...
Contact:

Plotting spectral profile

Post by Olivier Thizy »

Hello,


I have been working on a tool (ie: function) in MatLab to display spectral profile from FITS file we use to work with in astronomical spectroscopical community (ie: BeSS file format).

This is highly preliminary but just to start discussion around this subject, I publish here my initial code. Again, this is Work In Progress but maybe it can help you.

Just copy/paste this code into "otz_plot.m" file within matlab working directory.

Then go in matlab to a directory with at least one FITS file profile and type in the console for exemple:

otz_plot -h : for help
otz_plot 'filename.fits' : display your spectrum
otz_plot 'all' : display all FITS file spectra in your directory (could be messy if lot of them!)
otz_plot 'Ha' 'filename.fits' : to display Halpha from your spectrum (you can use 'Hb', 'Hg', 'D' for doublet, 'Ha+He' for both Halpha & HeI line...
otz_plot 6502 6603 0 4 'filename.fits' : here you specify the wavelength & Intensity limits
otz_plot 'filename.fits' -w 'Ha' : display in radial velocities Halpha
otz_plot 'filename.fits' -w 'Ha' -w 'Hb' -w 'Hg' : display in radial velocities Halpha, Hbeta & Hgamma
...

This is really WIP but do not hesitate to post here if you have comments or would like to see improvements in this code...

Code: Select all

function otz_plot(varargin)
%function otz_plot(filename)
% otz_plot.m
%
% Display and create a PDF of a graph of an echelle spectrum
%
% How to use:
% otz_plot [[L1 L2] Y1 Y2] filename1 -x param1 -x param2 param3 ...
% with L1 & L2 the wavelength (X axis) limits; if none, plot all spectrum
% with Y1 Y2 the Y-axis limits; if none, fit the graph
% use filename 'all' to use all .fits|.fit files in working directory
% 
% v1.0 / 20160816 (c) Olivier Thizy - published in ARAS forum
% Developped with MatLab R2016a
%

% TODO LIST for future versions !!!
% ---------------------------------
% lot of ideas to come up! :-)
%

% PARAMETERS
%-----------

% clean up some parameters
clear x;
clear y;
clear filename;
clear wavelengths;

% Option: -h or -help
% ==>print help text
param_h = 0; % 0: doesn't print help
%
% Option: -f filaname1 -f filename2 -f filename3...
% ==>allow multiple spectra on same graph
NbFiles = 1; % always one filename, first one doesn't need -f tag...
%
% Option: -c r|b|rainbow
% ==>color of the graph.
% r, b...: red, blue; based on MatLab colors
% rainbow: one color for each spectrum
param_c = 'b'; % plot spectrum in red
%
% Option: -p pdf|png [X Y]
% ==>print PDF, PNG or none graph with X/Y size if provided (not
% implemented yet!)
param_p = 'none'; % print PDF file
%
% Option: -m linetype
% ==>mark set of lines marks on the graph
% Np for planetary nebulae
% Be for Be stars
% CV for Cataclysmic Variables
% O for O stars, B for B stars, A for A stars, F, G, K, M
% WN for Wolf-Rayet WN type
% WC for Wolf-Rayet WC type
param_m = 'none';
%
% Option: -t title
% ==>type of title displayed: short, long (default), object, none
param_t = 'short';
%
% Option: -e #
% ==>display echelle spectrum in multiple graphs (#: Nb of graphs)
% ==>if #<0, display the full spectrum on top line then the -# graphs
%
% Option: -s f|t
% ==>allow shift for each graphs either fixed (f) or time based (t)
%
% Option: -y X1 X2
% ==>rescale each graph using mean of intensity (Y axis) within (X1,X2: wavelength or velocities) mean
param_scale = 0;
%
% Option: -w wavelength
% ==>plot in Radial Velocity using lambda wavelength
% ==> if used multiple times, plot different lines on same graph
c = double(300000); % speed of light
param_w = 0;
DefaultRV = 500; % graph will be displayed by default within +/-DefaultRV limits
%
% Option: -z
% ==>correct from heliocentric correction (radial velocity)
param_z = 0;
%

% Get parameters input
nVarargs = length(varargin);
Param = '';

Xlimit1 = -inf;
Xlimit2 = inf;
Ylimit1 = -inf;
Ylimit2 = inf;

NbFiles = 1;

% List of known lines (TODO: accuracy to be improved)
KnownLines = {'Ha','6563.0','50';'Hb','4861.0','50';'Hg','4340.0','50';'Hd','4101.0','50';
    'Ha+He','6610', '100';'H&K','3950','70';'HeI', '6678','50';'D','5888','38'};


% Get first arguments which do not require a command -X before
LastArg = nVarargs;
for j = nVarargs:-1:1
    if strcmp(varargin{j}(1),'-')
        LastArg = j-1;
    end
end

switch LastArg
    case 0
        return;
    case 1
        tmp = varargin{1};
        if  strcmp(tmp,'-h')|strcmp(tmp,'-help') % Special for HELP
            display('Tool to plot astronomical spectra in FITS & BeSS standard.');
            display('(c) Olivier Thizy / thizy@free.fr');
            display('Exemples:');
            display('>otz_plot ''filename.fits''');
            display('>otz_plot Ha ''filename.fits''');
            display('(here you can use: Ha, Hb, Hg, Hd, Ha+He, He, D, HK)');
            display('>otz_plot 6500 6600 ''filename.fits''');
            display('>otz_plot 6500 6600 0 5 ''filename.fits''');
            display('Sorry but extended HELP is not implemented yet...');
            return;
        else
            filename{1} = varargin{1};
        end
        NextVarN = 1;
    case 2
        Param = varargin{1};
        filename{1} = varargin{2};
        NextVarN = 2;
    case 3
        if varargin{2}(1) == '-'
            filename{1} = varargin{1};
            NextVarN = 1;
        else
            Xlimit1 = str2num(varargin{1});
            Xlimit2 = str2num(varargin{2});
            filename{1} = varargin{3};
            NextVarN = 3;
        end
    case 4
        Param = varargin{1};
        if varargin{3}(1) == '-'
            filename{1} = varargin{2};
            NextVarN = 2;
        else
            Ylimit1 = str2num(varargin{2});
            Ylimit2 = str2num(varargin{3});
            filename{1} = varargin{4};
            NextVarN = 4;
        end
    case 5
        Xlimit1 = str2num(varargin{1});
        Xlimit2 = str2num(varargin{2});
        Ylimit1 = str2num(varargin{3});
        Ylimit2 = str2num(varargin{4});
        filename{1} = varargin{5};
        NextVarN = 5;
    otherwise
        Xlimit1 = str2num(varargin{1});
        Xlimit2 = str2num(varargin{2});
        Ylimit1 = str2num(varargin{3});
        Ylimit2 = str2num(varargin{4});
        filename{1} = varargin{5};
        NextVarN = 6;
end

% get all FITS file from current directory if filename parameter is 'all'
if strcmp(filename{1},'all')
    location = '.\';
    ds = datastore(location,'Type','image','IncludeSubfolders',true,'FileExtensions',{'.fit','.fits','.FIT','.FITS'});
    % if multiple directories, can use {location1,location2,location3} syntax
    NbFiles = length(ds.Files);
    for k = 1:NbFiles
        filename{k} = char(ds.Files(k));
    end
end

% if first filename is '?', then open a dialog box to get it
if filename{1}=='?'
    [FileName, PathName] = uigetfile ('*.fits');
    filename{1} = [PathName, FileName];
end


% Parse the -X arguments
for k = NextVarN:nVarargs
    command=sscanf(varargin{k},'-%c');
    switch command
        case 'f' % multiple files plot
            NbFiles = NbFiles + 1;
            filename{NbFiles} = varargin{k+1};
            k=k+1;
        case 't' % title type
            param_t = varargin{k+1};
            k=k+1;
        case 'c' % color for the graph(s)
            param_c = varargin{k+1};
            k=k+1;
        case 'w' % Plot in velocity centered on wavelength(s)
            param_w = param_w + 1;
            if (k == nVarargs)
                if strcmp(Param,'')
                    wavelengths{param_w} = (Xlimit1 + Xlimit2)/2;                    
                else
                    for j = 1:size(KnownLines,1)
                        if strcmp(Param,KnownLines(j,1))
                            wavelengths{param_w} = str2double(KnownLines(j,2));
                        end
                    end            
                end
                k=k+1;
            else
                if strcmp(varargin{k+1}(1),'-')
                    wavelengths{param_w} = (Xlimit1 + Xlimit2)/2;
                    for j = 1:size(KnownLines,1)
                        if strcmp(Param,KnownLines(j,1))
                            wavelengths{param_w} = str2double(KnownLines(j,2));
                        end
                    end
                else
                    wavelengths{param_w} = str2double(varargin{k+1});
                    for j = 1:size(KnownLines,1)
                        if strcmp(varargin{k+1},KnownLines(j,1))
                            wavelengths{param_w} = str2double(KnownLines(j,2));
                            if Xlimit1 == -inf
                                Xlimit1 = double(-DefaultRV);
                            end
                            if Xlimit2 == inf
                                Xlimit2 = double(+DefaultRV);
                            end
                        end
                    end            
                end
                k=k+1;
            end
    end
end

% check if the line center is given in text/alias form
for k = 1:size(KnownLines,1)
    if strcmp(Param,KnownLines(k,1))
        if param_w == 0
            tmp1 = str2double(KnownLines(k,2));
            tmp2 = str2double(KnownLines(k,3));
            Xlimit1 = tmp1-tmp2;
            Xlimit2 = tmp1+tmp2;        
        else
            Xlimit1 = double(-DefaultRV);
            Xlimit2 = double(+DefaultRV);
        end
    end
end


% Browse all files selected
for k = 1:NbFiles
    % Extract FITS header and key infos from there
    FITSHeader = fitsinfo(filename{k});
    Data = fitsread(filename{k});
    % 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
        kwd = FITSHeader.PrimaryData.Keywords{l,1};
        
        switch kwd
            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
    
    % Create spectrum data: wavelength, intensity
    for a = 1 : NbPoints
        Lambda = Origine + double(a-1) * Pas;
        x (k,a) = double(Lambda);
        y (k,a) = Data (1,a);
    end
end

% display graph / open an half screen figure
figure('units','normalized','outerposition',[0.1 0.1 .5 .5])

% For each spectrum
for k = 1:NbFiles
    if param_w == 0
        plot(x(k,:), y(k,:), param_c);
        hold on;
    else
        % For each central wavelength
        for nw = 1:param_w
            for j = 1:size(x,2)
                w(k,j) = (x(k,j)-wavelengths{nw})*c/wavelengths{nw};
            end
            plot(w(k,:), y(k,:), param_c);
            hold on;
        end
    end
end

% Define X & Y limits for the graph
x_limit = [Xlimit1 Xlimit2];
xlim(x_limit)
y_limit = [Ylimit1 Ylimit2];
ylim(y_limit)

% Display title above the first graph with last file data
switch param_t
    case 'long'
        Titre2 = [Objet, ' (', Date, ' / ', ExpTime, ') @ ', SiteObs, ' w/ ', Instrument,'; (c) ', Observateur];
    case 'short'
        Titre2 = [Objet, ' - ', Observateur];
    case 'object'
        Titre2 = Objet;
    case 'none'
        Titre2 = '';
end

if strcmp(Titre2,'')
else
    title(Titre2, 'Interpreter','none');
end

if param_w == 0
    xlabel('Wavelenght in Â');
else
    xlabel('Velocity in km/s');
end
ylabel('Relative flux');


% save graph in PDF & PNG using 1st filename
[pathstr,name,ext] = fileparts(filename{1});
orient landscape;
outputfile = ['graph' name];
switch param_p
    case 'pdf'
        print (outputfile, '-painters', '-dpdf', '-bestfit');
    case 'png'
        orient portrait;
        print (outputfile, '-dpng');
end

% That's all folks !

Olivier Thizy
Posts: 370
Joined: Sat Sep 24, 2011 10:52 am
Location: in the french Alps...
Contact:

Re: Plotting spectral profile

Post by Olivier Thizy »

Hello,


here is v1.1 of the script with legend added:


Without -c option, the script is using standard colors for spectra in that order: b=blue, r=red, g=green, m=lagenta, y=yellow, c=cyan, k=black
You can specify colors with -c 'rgb' for exemple and graphs will be red, then green then blue and cycling...

I also added legends with DATE-OBS and central wavelength for plot in radial velocities.

Here is an exemple of three spectra of QR Vul I took recently, plotted in one line using command:
otz_plot 'Ha' 0.8 1.2 'all' -w
QR Vul
QR Vul
untitled.png (44.43 KiB) Viewed 8456 times

Code: Select all

function otz_plot(varargin)
%function otz_plot(filename)
% otz_plot.m
%
% Display and create a PDF of a graph of an echelle spectrum
%
% How to use:
% otz_plot [[L1 L2] Y1 Y2] filename1 -x param1 -x param2 param3 ...
% with L1 & L2 the wavelength (X axis) limits; if none, plot all spectrum
% with Y1 Y2 the Y-axis limits; if none, fit the graph
% use filename 'all' to use all .fits|.fit files in working directory
% 
% v1.0 / 20160816 (c) Olivier Thizy - published in ARAS forum
% v1.1 / 20160817 - modified color handling + legend added, published on ARAS
%
% Developped with MatLab R2016a
%

% TODO LIST for future versions !!!
% ---------------------------------
% lot of ideas to come up! :-)
%

% PARAMETERS
%-----------

% clean up some parameters
clear x;
clear y;
clear filename;
clear wavelengths;

% Option: -h or -help
% ==>print help text
param_h = 0; % 0: doesn't print help
%
% Option: -f filaname1 -f filename2 -f filename3...
% ==>allow multiple spectra on same graph
NbFiles = 1; % always one filename, first one doesn't need -f tag...
%
% Option: -c r|b|rainbow
% ==>color of the graph.
% r, b...: red, blue; based on MatLab colors
% rainbow: one color for each spectrum
param_c = 'brgmyck'; % plot spectra using this color scheme
% (b=blue, r=red, g=green, m=lagenta, y=yellow, c=cyan, k=black)
%
% Option: -p pdf|png [X Y]
% ==>print PDF, PNG or none graph with X/Y size if provided (not
% implemented yet!)
param_p = 'none'; % print PDF file
%
% Option: -m linetype
% ==>mark set of lines marks on the graph
% Np for planetary nebulae
% Be for Be stars
% CV for Cataclysmic Variables
% O for O stars, B for B stars, A for A stars, F, G, K, M
% WN for Wolf-Rayet WN type
% WC for Wolf-Rayet WC type
param_m = 'none';
%
% Option: -t title
% ==>type of title displayed: short, long (default), object, none
param_t = 'short';
%
% Option: -e #
% ==>display echelle spectrum in multiple graphs (#: Nb of graphs)
% ==>if #<0, display the full spectrum on top line then the -# graphs
%
% Option: -s f|t
% ==>allow shift for each graphs either fixed (f) or time based (t)
%
% Option: -y X1 X2
% ==>rescale each graph using mean of intensity (Y axis) within (X1,X2: wavelength or velocities) mean
param_scale = 0;
%
% Option: -w wavelength
% ==>plot in Radial Velocity using lambda wavelength
% ==> if used multiple times, plot different lines on same graph
c = double(300000); % speed of light
param_w = 0;
DefaultRV = 1000; % graph will be displayed by default within +/-DefaultRV limits
%
% Option: -z
% ==>correct from heliocentric correction (radial velocity)
param_z = 0;
%

% Get parameters input
nVarargs = length(varargin);
Param = '';

Xlimit1 = -inf;
Xlimit2 = inf;
Ylimit1 = -inf;
Ylimit2 = inf;

NbFiles = 1;

% List of known lines (TODO: accuracy to be improved)
KnownLines = {'Ha','6563.0','50';'Hb','4861.0','50';'Hg','4340.0','50';'Hd','4101.0','50';
    'Ha+He','6610', '100';'H&K','3950','70';'HeI', '6678','50';'D','5888','38'};


% Get first arguments which do not require a command -X before
LastArg = nVarargs;
for j = nVarargs:-1:1
    if strcmp(varargin{j}(1),'-')
        LastArg = j-1;
    end
end

switch LastArg
    case 0
        return;
    case 1
        tmp = varargin{1};
        if  strcmp(tmp,'-h')|strcmp(tmp,'-help') % Special for HELP
            display('Tool to plot astronomical spectra in FITS & BeSS standard.');
            display('(c) Olivier Thizy / thizy@free.fr');
            display('Exemples:');
            display('>otz_plot ''filename.fits''');
            display('>otz_plot Ha ''filename.fits''');
            display('(here you can use: Ha, Hb, Hg, Hd, Ha+He, He, D, HK)');
            display('>otz_plot 6500 6600 ''filename.fits''');
            display('>otz_plot 6500 6600 0 5 ''filename.fits''');
            display('Sorry but extended HELP is not implemented yet...');
            return;
        else
            filename{1} = varargin{1};
        end
        NextVarN = 1;
    case 2
        Param = varargin{1};
        filename{1} = varargin{2};
        NextVarN = 2;
    case 3
        if varargin{2}(1) == '-'
            filename{1} = varargin{1};
            NextVarN = 1;
        else
            Xlimit1 = str2num(varargin{1});
            Xlimit2 = str2num(varargin{2});
            filename{1} = varargin{3};
            NextVarN = 3;
        end
    case 4
        Param = varargin{1};
        if varargin{3}(1) == '-'
            filename{1} = varargin{2};
            NextVarN = 2;
        else
            Ylimit1 = str2num(varargin{2});
            Ylimit2 = str2num(varargin{3});
            filename{1} = varargin{4};
            NextVarN = 4;
        end
    case 5
        Xlimit1 = str2num(varargin{1});
        Xlimit2 = str2num(varargin{2});
        Ylimit1 = str2num(varargin{3});
        Ylimit2 = str2num(varargin{4});
        filename{1} = varargin{5};
        NextVarN = 5;
    otherwise
        Xlimit1 = str2num(varargin{1});
        Xlimit2 = str2num(varargin{2});
        Ylimit1 = str2num(varargin{3});
        Ylimit2 = str2num(varargin{4});
        filename{1} = varargin{5};
        NextVarN = 6;
end

% get all FITS file from current directory if filename parameter is 'all'
if strcmp(filename{1},'all')
    location = '.\';
    ds = datastore(location,'Type','image','IncludeSubfolders',false,'FileExtensions',{'.fit','.fits','.FIT','.FITS'});
%    ds = datastore(location,'Type','image','IncludeSubfolders',true,'FileExtensions',{'.fit','.fits','.FIT','.FITS'});
    % if multiple directories, can use {location1,location2,location3} syntax
    NbFiles = length(ds.Files);
    for k = 1:NbFiles
        filename{k} = char(ds.Files(k));
    end
end

% if first filename is '?', then open a dialog box to get it
if filename{1}=='?'
    [FileName, PathName] = uigetfile ('*.fits');
    filename{1} = [PathName, FileName];
end


% Parse the -X arguments
for k = NextVarN:nVarargs
    command=sscanf(varargin{k},'-%c');
    switch command
        case 'f' % multiple files plot
            NbFiles = NbFiles + 1;
            filename{NbFiles} = varargin{k+1};
            k=k+1;
        case 't' % title type
            param_t = varargin{k+1};
            k=k+1;
        case 'c' % color for the graph(s)
            param_c = varargin{k+1};
            k=k+1;
        case 'w' % Plot in velocity centered on wavelength(s)
            param_w = param_w + 1;
            if (k == nVarargs)
                if strcmp(Param,'')
                    wavelengths{param_w} = (Xlimit1 + Xlimit2)/2;                    
                else
                    for j = 1:size(KnownLines,1)
                        if strcmp(Param,KnownLines(j,1))
                            wavelengths{param_w} = str2double(KnownLines(j,2));
                        end
                    end            
                end
                k=k+1;
            else
                if strcmp(varargin{k+1}(1),'-')
                    wavelengths{param_w} = (Xlimit1 + Xlimit2)/2;
                    for j = 1:size(KnownLines,1)
                        if strcmp(Param,KnownLines(j,1))
                            wavelengths{param_w} = str2double(KnownLines(j,2));
                        end
                    end
                else
                    wavelengths{param_w} = str2double(varargin{k+1});
                    for j = 1:size(KnownLines,1)
                        if strcmp(varargin{k+1},KnownLines(j,1))
                            wavelengths{param_w} = str2double(KnownLines(j,2));
                            if Xlimit1 == -inf
                                Xlimit1 = double(-DefaultRV);
                            end
                            if Xlimit2 == inf
                                Xlimit2 = double(+DefaultRV);
                            end
                        end
                    end            
                end
                k=k+1;
            end
    end
end

% check if the line center is given in text/alias form
for k = 1:size(KnownLines,1)
    if strcmp(Param,KnownLines(k,1))
        if param_w == 0
            tmp1 = str2double(KnownLines(k,2));
            tmp2 = str2double(KnownLines(k,3));
            Xlimit1 = tmp1-tmp2;
            Xlimit2 = tmp1+tmp2;        
        else
            Xlimit1 = double(-DefaultRV);
            Xlimit2 = double(+DefaultRV);
        end
    end
end


% Browse all files selected
for k = 1:NbFiles
    % Extract FITS header and key infos from there
    FITSHeader = fitsinfo(filename{k});
    Data = fitsread(filename{k});
    % 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
        kwd = FITSHeader.PrimaryData.Keywords{l,1};
        
        switch kwd
            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
    
    dateGraph{k} = Date;
    
    % Create spectrum data: wavelength, intensity
    for a = 1 : NbPoints
        Lambda = Origine + double(a-1) * Pas;
        x (k,a) = double(Lambda);
        y (k,a) = Data (1,a);
    end
end

% display graph / open an half screen figure
figure('units','normalized','outerposition',[0.1 0.1 .5 .5])

NGraph = 1;

% For each spectrum
for k = 1:NbFiles
    if param_w == 0
        Graph{NGraph} = plot(x(k,:), y(k,:), param_c(1+rem(NGraph,size(param_c,2))));
        GraphLegend{NGraph} = dateGraph{k};
        NGraph = NGraph + 1;
        hold on;
    else
        % For each central wavelength
        for nw = 1:param_w
            for j = 1:size(x,2)
                w(k,j) = (x(k,j)-wavelengths{nw})*c/wavelengths{nw};
            end
            Graph{NGraph} = plot(w(k,:), y(k,:), param_c(1+rem(NGraph,size(param_c,2))));
            GraphLegend{NGraph} = [dateGraph{k} ' @ ' num2str(wavelengths{nw},'%d')];
            NGraph = NGraph + 1;
            hold on;
        end
    end
end

% Define X & Y limits for the graph
x_limit = [Xlimit1 Xlimit2];
xlim(x_limit)
y_limit = [Ylimit1 Ylimit2];
ylim(y_limit)

legend(GraphLegend);

% Display title above the first graph with last file data
switch param_t
    case 'long'
        Titre2 = [Objet, ' (', Date, ' / ', ExpTime, ') @ ', SiteObs, ' w/ ', Instrument,'; (c) ', Observateur];
    case 'short'
        Titre2 = [Objet, ' - ', Observateur];
    case 'object'
        Titre2 = Objet;
    case 'none'
        Titre2 = '';
end

if strcmp(Titre2,'')
else
    title(Titre2, 'Interpreter','none');
end

if param_w == 0
    xlabel('Wavelenght in Â');
else
    xlabel('Velocity in km/s');
end
ylabel('Relative flux');


% save graph in PDF & PNG using 1st filename
[pathstr,name,ext] = fileparts(filename{1});
orient landscape;
outputfile = ['graph' name];
switch param_p
    case 'pdf'
        print (outputfile, '-painters', '-dpdf', '-bestfit');
    case 'png'
        orient portrait;
        print (outputfile, '-dpng');
end

% That's all folks !

Cordialement,
Olivier Thizy
Olivier Thizy
Posts: 370
Joined: Sat Sep 24, 2011 10:52 am
Location: in the french Alps...
Contact:

otz_plot.m: tool to plot spectral profile(s)

Post by Olivier Thizy »

Hello,


here is my v1.5 code of my astronomical spectra plotting tool for MatLab. Several improvements have been done and I am using it all the time! :-)

Features:
-quick plotting, for exemple: otz_plot Ha 'file.fits'
-routine to open a FITS BeSS file and convert observer, location, object codes to more standard ones (ie: some cleaning!)...
(would certainly need to be improved depending on which spectra you are working with)
-multiple spectra plotting
-multiple wavelength plotting in velocity mode; ex. to quickly compare Halpha & Hbeta: otz_plot 'file.fits' -w Ha -w Hb
-intensity rescaling & time series, periodic or not (next step would certainly be surface plots!)
-control colors, line styles, legend(s), title with quick codes
-export in PNG or PDF format if needed
-use direct link to CDS to get data of the object(s) such as Radial Velocities, RA/Dec
-correct from object Radial Velocity (catalog) and Heliocentric RV (calculated)

It is fairly basic but I have written for my own purpose to quickly do nice graphics with spectra...


Installation:
You need four files which I all list in this message below, some functions can be used independantly of course... Those files are:
otz_plot.m : the main function to plot spectra
rfitsinfo.m : get header & data from an astronomical spectrum (single HDU) FITS file
CDS_Infos.m : get data from CDS server
otz_hjd.m : calculate heliocentric julian date & RV correction to apply

Copy the files in your MatLab directory (Documents\MATLAB for Windows)
and type otz_plot -h to get this help file which provides several examples:
Tool to plot astronomical spectra in FITS & BeSS standard.
(c) Olivier Thizy / thizy@free.fr
v1.5 - 20150909 - tested with MatLab R2016a
---------------------------------
>otz_plot -h
Will give you this help...
>otz_plot -lines
Display the list of known lines with abreviations
----- Options to select the files to display -----
>otz_plot 'filename.fits'
>otz_plot all
>otz_plot 'asdb_vvcep*.fits'
>otz_plot ?
>otz_plot 'filename.fits' -f 'file2.fits' -f 'file3.fits'
----- Options to select the spectral domain to display -----
>otz_plot 6500 6600 'filename.fits'
>otz_plot 6500 6600 0 5 'filename.fits'
(here 0 5 are the Y-axis limits for the graphs)
>otz_plot Ha 'filename.fits'
(some codes can be used to quickly select the line to display; 'otz_plot -l' to list the line codes)
>otz_plot 'filename.fits' -w Ha -w Hb -w Hg
(the same line codes can be used for velocity plots)
-w (wavelength) options are the same as the line codes given by 'otz_plot -l'
----- Options for line color(s) & style(s) -----
>otz_plot 6500 6600 0 5 'all' -c bgr
-c (color) options: b=blue, r=red, g=green, m=lagenta, y=yellow, c=cyan, k=black
>otz_plot 6500 6600 0 5 'all' -d sp
-d (line style) options: s=solid, d=dash, p=dotted, m=dash-dot
----- Options for legend(s) -----
>otz_plot 'filename.fits' -l X
-l (legend): you can specify what data you want to display in legend for each graph
'legend' are codes also used for the -t (title) option; they are:
X or x means no legend
D: date in DD-MMM-YYYY HH:MM:SS
e: exposure in sec)
i: instrument
h: heliocentric JD -2400000
j: julian date -2400000
J: full julian date + 3 digits
l: location site code
L: full location name
n: filename
o: Observer code; O: Observer full name
r: resolution power of the instrument
t: target name (standard)
T: object name from FITS file
w: wavelength, no digit (when option -w is used)
Exemple: '-l Do' to display both date & observer code in legend
(by default legend option is 'Dwo': Date - wavelength - observer code
----- Options for main title -----
>otz_plot 'filename.fits' -t o
-t (title): you can specify what data you want to display as title
X or x means no legend; 'title' is similar to -l (legend) options
(by default title option is 'tDeLiO': target / Date / exposure time / full location / instrument / full observer name
----- Options scaling Y-axis & X-axis -----
otz_plot -500 600 0 10 all -s 6540 6542 -w Ha
-s: rescale to 1 beteen 6540A & 6542A
>otz_plot 'filename.fits' -x -18.3 -f 'filename2.fits' -x 23
-x RV: use the object RV if known (see rfitsinfo.m / use CDS to get RV)
-x H: heliocentric correction (minor bug in calculation)
-x RVH: object RV *and* heliocentric correction
-x nA: apply a shift of 'n' Angstroems
-x v: apply a shift of 'v' in km/s
otz_plot -500 600 0 10 all -t object -y -0.01 -w Ha
-y i|t|p scale (epoch period) with options:
-y i scale / intensity plot; scale = intensity shift between each graph (ex: 0.5)
-y t scale / time plot; scale = days to intensity scaling factor (0.01 for exemple; 100days = 1 in intensity)
-y p scale epoch period / phase plot; phase to intensity scaling factor, epoch in Heliocentric Julian Day, period in days
(for ex VV Cep: -y p 100 2438481 7430) [=Emin / Wright, R.A.S.C. Jour., Vol 71, 1977]
(note that -y option is usually used in conjunction with -s option to ensure proper intensity scale)
otz_plot Ha all -t object -s 6540 6542 -y 0.5 -w Ha
-y: if positive, shifts all spectra by 0.5 Intensity
---------------------------------
I hope you'll enjoy it !


The four functions are in the ZIP file:
otz_plot_v1-5_20160909.zip
(17.65 KiB) Downloaded 645 times

Cordialement,
Olivier Thizy
Post Reply