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 !
