Page 1 of 1

### Plotting spectral profile Posted: Tue Aug 16, 2016 9:02 pm
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 parametersclear x;clear y;clear filename;clear wavelengths;% Option: -h or -help% ==>print help textparam_h = 0; % 0: doesn't print help%% Option: -f filaname1 -f filename2 -f filename3...% ==>allow multiple spectra on same graphNbFiles = 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 spectrumparam_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 typeparam_m = 'none';%% Option: -t title% ==>type of title displayed: short, long (default), object, noneparam_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) meanparam_scale = 0;%% Option: -w wavelength% ==>plot in Radial Velocity using lambda wavelength% ==> if used multiple times, plot different lines on same graphc = double(300000); % speed of lightparam_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 inputnVarargs = 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 beforeLastArg = nVarargs;for j = nVarargs:-1:1    if strcmp(varargin{j}(1),'-')        LastArg = j-1;    endendswitch 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));    endend% if first filename is '?', then open a dialog box to get itif filename{1}=='?'    [FileName, PathName] = uigetfile ('*.fits');    filename{1} = [PathName, FileName];end% Parse the -X argumentsfor 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    endend% check if the line center is given in text/alias formfor 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    endend% Browse all files selectedfor 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);    endend% display graph / open an half screen figurefigure('units','normalized','outerposition',[0.1 0.1 .5 .5])% For each spectrumfor 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    endend% Define X & Y limits for the graphx_limit = [Xlimit1 Xlimit2];xlim(x_limit)y_limit = [Ylimit1 Ylimit2];ylim(y_limit)% Display title above the first graph with last file dataswitch 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 = '';endif strcmp(Titre2,'')else    title(Titre2, 'Interpreter','none');endif param_w == 0    xlabel('Wavelenght in Â');else    xlabel('Velocity in km/s');endylabel('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 !`

### Re: Plotting spectral profile Posted: Thu Aug 18, 2016 5:14 pm
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
untitled.png (44.43 KiB) Viewed 4536 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 parametersclear x;clear y;clear filename;clear wavelengths;% Option: -h or -help% ==>print help textparam_h = 0; % 0: doesn't print help%% Option: -f filaname1 -f filename2 -f filename3...% ==>allow multiple spectra on same graphNbFiles = 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 spectrumparam_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 typeparam_m = 'none';%% Option: -t title% ==>type of title displayed: short, long (default), object, noneparam_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) meanparam_scale = 0;%% Option: -w wavelength% ==>plot in Radial Velocity using lambda wavelength% ==> if used multiple times, plot different lines on same graphc = double(300000); % speed of lightparam_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 inputnVarargs = 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 beforeLastArg = nVarargs;for j = nVarargs:-1:1    if strcmp(varargin{j}(1),'-')        LastArg = j-1;    endendswitch 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));    endend% if first filename is '?', then open a dialog box to get itif filename{1}=='?'    [FileName, PathName] = uigetfile ('*.fits');    filename{1} = [PathName, FileName];end% Parse the -X argumentsfor 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    endend% check if the line center is given in text/alias formfor 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    endend% Browse all files selectedfor 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);    endend% display graph / open an half screen figurefigure('units','normalized','outerposition',[0.1 0.1 .5 .5])NGraph = 1;% For each spectrumfor 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    endend% Define X & Y limits for the graphx_limit = [Xlimit1 Xlimit2];xlim(x_limit)y_limit = [Ylimit1 Ylimit2];ylim(y_limit)legend(GraphLegend);% Display title above the first graph with last file dataswitch 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 = '';endif strcmp(Titre2,'')else    title(Titre2, 'Interpreter','none');endif param_w == 0    xlabel('Wavelenght in Â');else    xlabel('Velocity in km/s');endylabel('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

### otz_plot.m: tool to plot spectral profile(s) Posted: Fri Sep 09, 2016 10:13 pm
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