% PLOTTSVSC   Tool to plot variable files of the TSVSC catalog.
%    PLOTTSVSC(filename)
%
% example plottsvsc('TN-N222212012-9-67-2.txt')

function plottsvsc(filename)

% === Read TAROT variable file
fclose('all');
[fid msg]=fopen(filename,'r');
if (fid==-1)
	disp(msg);
end
ndata=-1;
while (feof(fid)==0)
	ligne=fgets(fid);
	if (ligne(1:3)=='END')
		ndata=0;
		continue;
	end
	if (ndata==-1)
		ligne=ligne(1:length(ligne)-2);
		[keyword ligne]=strtok(ligne,'=');
		keyword=strtok(keyword);
		[dummy ligne]=strtok(ligne);
		value=regexprep(ligne,' ','');
		if (strcmp(keyword,'NAME')==1)
			name=value;
		end
		if (strcmp(keyword,'RA')==1)
			ra=str2num(value);
		end
		if (strcmp(keyword,'DEC')==1)
			dec=str2num(value);
		end
		if (strcmp(keyword,'FILTER')==1)
			filter=value;
		end
		if (strcmp(keyword,'PERIOD')==1)
			period=str2num(value);
		end
		if (strcmp(keyword,'JDPHASE')==1)
			jdphase=str2num(value);
		end
		if (strcmp(keyword,'C1')==1)
			c1=str2num(value);
		end
		for k=1:5
			if (strcmp(keyword,['A',num2str(k)])==1)
				a(k)=str2num(value);
			end
			if (strcmp(keyword,['B',num2str(k)])==1)
				b(k)=str2num(value);
			end
		end
	elseif (ndata>=0)
		dummy=sscanf(ligne,'%lf');
		ndata=ndata+1;
		jds(ndata)=dummy(1);
		mags(ndata)=dummy(2);
		dmags(ndata)=dummy(3);
		flags(ndata)=dummy(4);
	end
end
fclose(fid);

% === Conversion UT -> barycentric (low accuracy: about +/-60s)
dr=pi/180;
oneyear=365.2425;
jd0=2453636.0;
dut2tt=64.184/86400;
epsi=23.438689*dr;
UA=1.49597870e11;
CLIGHT=2.99792458e8;
xearth=cos(2*pi*(jds+dut2tt-jd0)/oneyear);
yearth=sin(2*pi*(jds+dut2tt-jd0)/oneyear);
zearth=0*xearth;
x=xearth;
y=yearth*cos(epsi)-zearth*sin(epsi);
z=yearth*sin(epsi)+zearth*cos(epsi);
xx=cos(ra*dr)*cos(dec*dr);
yy=sin(ra*dr)*cos(dec*dr);
zz=sin(dec*dr);
d=x.*xx+y.*yy+z.*zz;
djd=d*UA/CLIGHT/86400.;
jdbarys=jds+djd;

% === compute the fit phases
phasefits=linspace(0,1,100);
n=length(a);
magfits=c1;
for k=1:n
	magfits=magfits+a(k).*cos(2*pi*k*phasefits)+b(k).*sin(2*pi*k*phasefits);
end

% === compute the phases for observation dates
phases=(jdbarys-jdphase)/period;
phases=phases-floor(phases);

% === Plot mag=f(jd)
clf;
set(gcf,'PaperUnits','centimeters');
set(gcf,'Units','centimeters');
set(gcf,'PaperPosition',[0 0 10 5]); % position sur la feuille
set(gcf,'Position',[0 2 10 10]); % position sur l'ecran
set(gcf,'InvertHardcopy','off');
set(gcf,'color',[1 1 1]);
h=subplot('position',[0.15 0.12 0.8 0.8]);
plot(jdbarys-jdphase,mags,'.k');
hold on
plot([jdbarys-jdphase ; jdbarys-jdphase],[mags+dmags ; mags-dmags],'k');
xlabel(['Days from JD_{bary} ',num2str(jdphase)],'FontSize',10,'FontName','Times');
ylabel([filter,' magnitude'],'FontSize',10,'FontName','Times');
title([name,' (P=',num2str(period),' days)'],'FontSize',10,'FontName','Times');
set(gca,'Ydir','reverse','box','on','FontSize',10,'FontName','Times');
pause;

% === Plot mag=f(phase)
clf;
h=subplot('position',[0.15 0.12 0.8 0.8]);
plot(phases,mags,'.k');
hold on
plot([phases ; phases],[mags+dmags ; mags-dmags],'k');
plot(phasefits,magfits,'r','linewidth',2);
xlabel(['Phase (0 for JD_{bary} ',num2str(jdphase),')'],'FontSize',10,'FontName','Times');
set(gca,'Ydir','reverse','box','on','FontSize',10,'FontName','Times');
ylabel([filter,' magnitude'],'FontSize',10,'FontName','Times');
title([name,' (P=',num2str(period),' days)'],'FontSize',10,'FontName','Times');
