I wrote an ATEL about it - thanks to Steve Shore for helping me out:
http://www.astronomerstelegram.org/?read=1196
Here are my own spectra (all in BeSS database) and I encourage everyone to follow this target, specially during such outburst.
Using MatLab 3D script, I was able to do some 3D graphs as well as spectrogram (here they are spectra with the first one substracted, to better show the emission).
The emission is very bright on Halpha, well visible on Hbeta and can be detected on Hgamma. The He I 5876 and line at 4820A do not show emission but clear periodic behaviour (stellar rotation) and at the time of emission I seem to detect a change in the spectrum too.
Code for the MatLab Graph3D.m
%
% Graph2D
%
% Plot a graph from a serie of spectra
% X axis: wavelength (in radial velocity)
% Y axis: time, interpolated, in JD
%
% TODO
% ****
% 1/ fold date into phase using Period & HJD0
% 2/ move from JD to HJD
% 3/ correct wavelength from heliocentric RV
% 4/ dsplay dual labels on date axis: HJD-2457000 *and* Year
%
clear
% Create list of FITS files in the directpry & subdirectories
ds = datastore('.\','Type','image','IncludeSubfolders',false,'FileExtensions',{'.fit','.fits','.FIT','.FITS'});
% Global variables
c=300000.0; % light speed in km/s
JD0=2457000.0;
% Key data from target (V442 And) from CDS & litterature
Name='V442 And';
Author='(c) O.Thizy';
RATxt = '01 03 53.3583';
RA = 0.2788;
DecTxt = '+47 38 32.260';
Dec = 0.8315;
SpectralType = 'Be star';
MagB = 6.7900;
MagV = 6.8200;
MagR = 6.8200;
RV = -55.00;
HJD0=2457987.85;
Period=2.61507;
% Select the wavelength to display
Lambda0=4920.0; % Line near Hbeta
%=6678.0; % HeI singlet
%=6563.8; % Halpha
%=6678.0; % HeI singlet
%=4920.0; % Line near Hbeta
%=4861.0; % Hbeta
%=4340.0; % Hgamma
%=5876.0; % HeI triplet
%=4921.0; %HeI singlet
% Define which is the reference spectrum (usually the first one)
REFNum=1;
Titre=[Name ' / ' num2str(Lambda0) 'A / ' Author];
% Define the spectral domain to display
v1=-500.0; % in km/s
v2=+500.0; % in km/s
% Define the beginning of the spectral domain for all your spectra
% and the dispersion
%CRVAL=6507.5; % Valid for my echelle spectra only - Halpha!!!
CRVAL=4825.5; % Valid for my echelle spectra only - Hbeta!!!
%CRVAL=4277.5; % Valid for my echelle spectra only - Hgamma!!!
%CRVAL=5831.5; % Valid for my echelle spectra only - Na doublet!!!
CDELT=0.05; % idem
% Calculate the index position of the spectral domain within the
% interpolated spectrum
x1 = int32((Lambda0 + v1*Lambda0/c - CRVAL)/CDELT);
x2 = int32((Lambda0 + v2*Lambda0/c - CRVAL)/CDELT);
dX=int32(x2-x1)+1; % Spectrum now will go through (1:dX)
nX = int32(1000);
nY = 100;
REF = fitsread(ds.Files{REFNum}); % Reference spectrum vector (for substrated graph)
% NAXIS=63281; % Nb of Elements
MATRIX=zeros(length(ds.Files),dX);
MATRIX_SUB=zeros(length(ds.Files),dX);
MATRIX_PHASE=zeros(length(ds.Files),dX);
%OBS_DATE=zeros(length(ds.Files));
for i = 1:length(ds.Files)
fName = char(ds.Files(i));
H=fitsinfo(ds.Files{i});
Data = fitsread(ds.Files{i}); % spectrum vector
% searching for key header values
if max(size(H.Contents)) == 1
kwd = H.PrimaryData.Keywords ;
elseif max(size(H.Contents)) > 3
kwd = [H.PrimaryData.Keywords;H.Image(max(size(H.Contents))-2).Keywords];
else
kwd = H.PrimaryData.Keywords ;
end
HeaderMax = max(size(kwd));
for L=1:HeaderMax
k = kwd {L,1};
kc = kwd {L,2};
kt = kwd {L,3};
switch k
% CALIBRATION & REDUCTION data
case 'NAXIS1'
NAXIS1 = kwd {L,2};
case 'CRVAL1'
CRVAL1 = kwd {L,2};
case 'CDELT1'
CDELT1 = kwd {L,2};
case 'JD-MID'
JD = kwd {L,2};
case 'EXPTIME2' % Total exposure duration (including in-between read-out time)
ExpTimeTotal = kwd {L,2};
case 'DATE-OBS'
if length(kwd {L,2}) == 17
ObsDate = datetime(kwd {L,2},'InputFormat','uuuuMMdd''T''HH:mm:ss','TimeZone','UTC');
elseif length(kwd {L,2}) == 18
ObsDate = datetime(kwd {L,2},'InputFormat','uuuuMMdd''T:''HH:mm:ss','TimeZone','UTC');
elseif length(kwd {L,2}) == 19
if strcmp(kwd {L,2}(14),'-')
ObsDate = datetime(kwd {L,2},'InputFormat','uuuu-MM-dd''T''HH-mm:ss','TimeZone','UTC');
else
ObsDate = datetime(kwd {L,2},'InputFormat','uuuu-MM-dd''T''HH:mm:ss','TimeZone','UTC');
end
elseif length(kwd {L,2}) == 20
ObsDate = datetime(kwd {L,2},'InputFormat','uuuu-MM-dd''T''HH:mm:ss.','TimeZone','UTC');
elseif length(kwd {L,2}) == 21
ObsDate = datetime(kwd {L,2},'InputFormat','uuuu-MM-dd''T''HH:mm:ss.S','TimeZone','UTC');
elseif length(kwd {L,2}) == 22
ObsDate = datetime(kwd {L,2},'InputFormat','uuuu-MM-dd''T''HH:mm:ss.SS','TimeZone','UTC');
elseif length(kwd {L,2}) == 30
ObsDate = datetime(kwd {L,2},'InputFormat','uuuu-MM-dd''T''HH:mm:ss.SSSSSSSSSS','TimeZone','UTC');
elseif length(kwd {L,2}) == 31
ObsDate = datetime(kwd {L,2},'InputFormat','uuuu-MM-dd''T''HH:mm:ss.SSSSSSSSSSS','TimeZone','UTC');
elseif length(kwd {L,2}) == 32
ObsDate = datetime(kwd {L,2},'InputFormat','uuuu-MM-dd''T''HH:mm:ss.SSSSSSSSSSSS','TimeZone','UTC');
elseif length(kwd {L,2}) >= 23
ObsDate = datetime(kwd {L,2},'InputFormat','uuuu-MM-dd''T''HH:mm:ss.SSS','TimeZone','UTC');
elseif isempty(kwd {L,2})
ObsDate = kwd {L,2};
end
end
end
% Lambda = CRVAL1 + CDELT1 * x
% x = (Lambda - CRVAL1)/CDELT1
% v = (Lambda-Lambda0)/Lambda0*c
% Lambda = Lambda0 + v*Lambda0/c
%Xx1 = int32((Lambda0 + v1*Lambda0/c - CRVAL1)/CDELT1);
%Xx2 = int32((Lambda0 + v2*Lambda0/c - CRVAL1)/CDELT1);
% dXx = int32(Xx2-Xx1)+1;
% Xnew=linspace(1,dXx,nX);
OBS_DATE(i)=(JD-JD0);
% CalcJD = juliandate(ObsDate + seconds(ExpTimeTotal / 2))
[HJD,ObjVel]=otz_hjd(JD,[RA Dec]);
OBS_PHASE(i)=(HJD-JD0)/Period-fix((HJD-JD0)/Period);
% MATRIX(i,1:dX)=Data(x1:x2);
MATRIX_SUB(i,1:dX)=Data(x1:x2)/mean(Data(x2-50:x2-40))-REF(x1:x2)/mean(REF(x2-50:x2-40));
MATRIX(i,1:dX)=Data(x1:x2)/mean(Data(x2-50:x2-40));
end % end of loop for each file/spectrum
% pcolor(MATRIX), shading interp;
T=1:length(ds.Files);
Tnew=linspace(1,length(ds.Files),nY);
%X=double(1:dX);
%Vnew=linspace(double(1),double(dX),double(nX));
Lambda=CRVAL1+double((x1:x2))*CDELT1;
Velocity=(Lambda-Lambda0)/Lambda0*c;
Y=linspace(OBS_DATE(1),OBS_DATE(length(ds.Files)),nY);
Y_PHASE=linspace(0,1,nY);
IMAGE=interp1(T,MATRIX,Tnew);
IMAGE_SUB=interp1(T,MATRIX_SUB,Tnew);
%IMAGE2=interp2(T,X,MATRIX(:,x1:x2),Tnew,Vnew);
T=1:length(ds.Files);
Tnew=linspace(0,1,nY);
IMAGE_PHASE=interp1(T,MATRIX,Tnew);
% At the end, create & display:
%
% 2D image in Velocity Vs JD-JD0
% 2D image in Velocity Vs Phase [0:1]
% 2D image + 3D graph (spectra) in Velocity Vs JD-JD0
% 2D image + 3D graph (spectra) in Velocity Vs Phase [0:1]
% Same with REF spectra substracted
%
%pcolor(Velocity,Y,IMAGE),shading interp,colormap gray;
%pcolor(IMAGE),shading interp,colormap gray;
%figure;
%surf(Velocity,Y,IMAGE),shading interp,colormap parula,colorbar,ylabel('JD-2458000'),xlabel('velocity (km/s)');
% colormap jet, colormap gray, colormap parula
% Use plot for 1D, pcolor for 2D and surf for 3D, contourf for art !
%contourf(Velocity,Y,IMAGE),shading interp,colormap gray,colorbar,ylabel('JD-2458000'),xlabel('velocity (km/s)');
% Note: whos VARIABLE gives structure information on the variable !
% 2D image in Velocity Vs JD-JD0
figure
pcolor(Velocity,Y,IMAGE),shading interp,colormap gray,colorbar,ylabel(['JD-' num2str(JD0)]),xlabel('velocity (km/s)'),title(Titre);
% 2D image in Velocity Vs JD-JD0 with REF spectrum substracted
figure
pcolor(Velocity,Y,IMAGE_SUB),shading interp,colormap gray,colorbar,ylabel(['JD-' num2str(JD0)]),xlabel('velocity (km/s)'),title(Titre);
% Graph 2D with all spectra (surface plot & 3D plot, time shifted
figure;
hold on
pcolor(Velocity,Y,IMAGE),shading interp,colormap parula,colorbar,ylabel(['JD-' num2str(JD0)]),xlabel('velocity (km/s)'),title(Titre);
for i=1:length(OBS_DATE)
plot3(Velocity,ones(size(Velocity))*OBS_DATE(i),MATRIX(i,:))
end
% Graph 2D with all spectra (surface plot & 3D plot, folded in phased time
%figure
% hold on
% pcolor(Velocity,Y_PHASE,IMAGE_PHASE),shading interp,colormap parula,colorbar,ylabel('Phase'),xlabel('velocity (km/s)'),title(Titre);
%for i=1:length(OBS_DATE)
% plot3(Velocity,ones(size(Velocity))*OBS_PHASE(i),MATRIX(i,:))
%end
Aucun commentaire:
Enregistrer un commentaire