clear
clf

pixpermin=0:0.01:1;

pix=0.6054; %arcmin
rs=6.96e8;  %radius of sun in metres
pixm=pix*7.22e5; % 1pix in metres
v=pixpermin*pixm/60;
nHz=1e9*(v/(2*pi*rs)); % rotation freq in nHz
figure(1)
fh1=subplot(3,2,1);
set(fh1,'visible','off');
plot(pixpermin,v); grid;
title('velocity at photosphere: m.s^{-1} vs. pix.min^{-1} ');
xlabel('v/(pix/min)');
ylabel('v/(m.s^{-1})');

fh2=subplot(3,2,2);
set(fh2,'visible','off');
%a2=axes('parent',gcf,'YaxisLocation','right','position',get(fh2,'position'))
%a1=axes('parent',gcf,'YaxisLocation','left','position',get(fh2,'position'))
a2=axes('YaxisLocation','right','position',get(fh2,'position'))
a1=axes('YaxisLocation','left','position',get(fh2,'position'))


%plotyy(pixpermin,nHz,pixpermin,nHz); grid;
plot(a1,pixpermin,nHz); grid;
set(a2,'ylim',get(a1,'ylim'))
ylabels=str2num(get(a1,'yticklabel'));
set(a2,'yticklabel',num2str(1./(24*3600*(ylabels*1e-9)),2));
set(get(a2,'ylabel'),'string','period/days')
%add a second axis
title('velocity at photosphere: nHz at equator vs. pix.min^{-1} ');
xlabel('v/(pix/min)');
ylabel('rotation rate /nHz');



% do the same at higher resolution
pixpermin=0.181:0.001:0.3;
pixm=pix*7.22e5;
v=pixpermin*pixm/60;

fh3=subplot(3,2,3);
set(fh3,'visible','off');
plot(pixpermin,v); grid;
title('velocity at photosphere: m.s^{-1} vs. pix.min^{-1} ');
xlabel('v/(pix/min)');
ylabel('v/(m.s^{-1})');

fh4=subplot(3,2,4);
set(fh4,'visible','off');
a2=axes('YaxisLocation','right','position',get(fh4,'position'))
a1=axes('YaxisLocation','left','position',get(fh4,'position'))
set(a2,'xtick',[]);

nHz=1e9*(v/(2*pi*rs));
%plotyy(pixpermin,nHz,pixpermin,nHz); grid;
plot(a1,pixpermin,nHz); grid;
set(a2,'ylim',get(a1,'ylim'))
ylabels=str2num(get(a1,'yticklabel'));
set(a2,'yticklabel',num2str(1./(24*3600*(ylabels*1e-9)),2));
set(get(a2,'ylabel'),'string','period/days')
%add a second axis
title('velocity at photosphere: nHz at equator vs. pix.min^{-1} ');
xlabel('v/(pix/min)');
ylabel('rotation rate /nHz');

% ********************************************************************
% do the same at higher resolution
pixpermin=0:0.001:0.10;
pixm=pix*7.22e5;
v=pixpermin*pixm/60;

fh5=subplot(3,2,5);
set(fh5,'visible','off');
plot(pixpermin,v); grid;
title('velocity at photosphere: m.s^{-1} vs. pix.min^{-1} ');
xlabel('v/(pix/min)');
ylabel('v/(m.s^{-1})');

fh6=subplot(3,2,6);
set(fh6,'visible','off');
a2=axes('YaxisLocation','right','position',get(fh6,'position'))
a1=axes('YaxisLocation','left','position',get(fh6,'position'))
set(a2,'xtick',[]);
nHz=1e9*(v/(2*pi*rs));
%plotyy(pixpermin,nHz,pixpermin,nHz); grid;
plot(a1,pixpermin,nHz); grid;
set(a2,'ylim',get(a1,'ylim'))
ylabels=str2num(get(a1,'yticklabel'));
set(a2,'yticklabel',num2str(1./(24*3600*(ylabels*1e-9)),2));
set(get(a2,'ylabel'),'string','period/days')
%add a second axis
title('velocity at photosphere: nHz at equator vs. pix.min^{-1} ');
xlabel('v/(pix/min)');
ylabel('rotation rate /nHz');


nHz=1e9*(v/(2*pi*rs));
%plotyy(pixpermin,nHz,pixpermin,nHz); grid;
plot(a1,pixpermin,nHz); grid;
set(a2,'ylim',get(a1,'ylim'))
ylabels=str2num(get(a1,'yticklabel'));
set(a2,'yticklabel',num2str(1./(24*3600*(ylabels*1e-9)),3));
set(get(a2,'ylabel'),'string','period/days')
%add a second axis
title('velocity at photosphere: nHz at equator vs. pix.min^{-1} ');
xlabel('v/(pix/min)');
ylabel('rotation rate /nHz');