%---------------------------------------------------------------------------
%
%     MATLAB program to plot results of buckling.f 
%     See readme.txt and buckling.f for a description
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%
%   (c)2004,2012,2015 Maurice H.P.M. van Putten
%	This program is distributed under the GNU General Public Licence
%	You should have received a copy of the GNU General Public License
%	along with this program. If not, see <http://www.gnu.org/licenses/>.
%
%----------------------------------------------------------------------------

%Load output buckling.f
load fort.301;f31=fort;load fort.302;f32=fort;
load fort.303;f33=fort;load fort.304;f34=fort;load fort.305;f35=fort;
load fort.401;f41=fort;load fort.402;f42=fort;
load fort.403;f43=fort;load fort.404;f44=fort;load fort.405;f45=fort;

%-------------------------------------------------------------
% Fig. 10a, van Putten, 2012, Prog. Theor. Phys.,127, 331 
% Fig. 2, van Putten, 2002, ApJ, 575, L71
%-------------------------------------------------------------
figure(1);hold on;box on;%grid on;
plot(f31(:,1),f31(:,2));plot(f32(:,1),f32(:,2));
plot(f33(:,1),f33(:,2));plot(f34(:,1),f34(:,2));plot(f35(:,1),f35(:,2));
FS=18;
xlabel('\itb/R','FontSize',FS);ylabel('q_c','FontSize',FS)
x0=0.14;y0=0.03;
text(0.14,2.40,'5','FontSize',FS);text(0.19,2.38,'4','FontSize',FS)
text(0.39-x0,2.3+y0,'3','FontSize',FS);text(0.49-x0,2.15+y0,'2','FontSize',FS)
text(0.55-x0,1.875+y0,'m=1','FontSize',FS);axis([0 1 1.5 2.5])
delta=f31(:,1);
I=1:20:900;
T10=1;q=1.5+0.15./(delta/0.2).^2*T10;plot(delta(I),q(I),'--');
T10=2;q=1.5+0.15./(delta/0.2).^2*T10;plot(delta(I),q(I),'--');
T10=3;q=1.5+0.15./(delta/0.2).^2*T10;plot(delta(I),q(I),'--');
text(0.33,1.70,'T_{10}=1,2,3','FontSize',FS)

title('Stability diagram','FontSize',FS)
text(0.6,1.675,'Fig. 10a in van Putten, 2012','FontSize',12)  
text(0.6,1.625,'Prog. Theor. Phys., 127, 331','FontSize',12)  

print -deps2 pA.eps; print -djpeg pA.jpg;

%-------------------------------------------------------------
% Fig. 10b, van Putten, 2012, Prog. Theor. Phys.,127, 331 
%-------------------------------------------------------------
figure(2);
for P=1:5
    if(P==1)delta=f31;end;if(P==2)delta=f32;end
    if(P==3)delta=f33;end;if(P==4)delta=f34;end
    if(P==5)delta=f35;end
    j=0;
    for T10=0.1:0.1:4
        j=j+1;
        q=1.5+0.15./(delta(:,1)/0.2).^2*T10;
        k=0;
        for i=1:length(delta)
            if(delta(i,2)>q(i) & k==0)k=i;end;
        end
        t(j)=T10;
        dc(j)=delta(k);
        qc(j)=q(k);
    end
    plot(dc,t);hold on;

end
xlabel('\itb/R','FontSize',FS);
ylabel('Temperature [10^{10}K]','FontSize',FS);
text(0.200,1.5,'m=1','FontSize',FS);
text(0.200,2.1,'2','FontSize',FS);text(0.190,2.6,'3','FontSize',FS);
text(0.165,3.0,'4','FontSize',FS);text(0.145,3.3,'5','FontSize',FS);
axis([0 0.3 0 4]);%grid on;

title('Critical temperatures','FontSize',FS)
text(0.19,0.4,'Fig. 10b in van Putten, 2012','FontSize',12)  
text(0.19,0.2,'Prog. Theor. Phys., 127, 331','FontSize',12)  
 
print -depsc pB.eps;print -djpeg pB.jpg

%-----------------------------------------------------------
% Fig. 3, van Putten, 2002, ApJ, 575, L71
%-------------------------------------------------------------
figure(3);hold on;%grid on;
I=1:20:900;
plot(f41(:,1),f41(:,5));plot(f42(:,1),f42(:,5));
plot(f43(:,1),f43(:,5));plot(f44(:,1),f44(:,5));
plot(f45(:,1),f45(:,5));
plot(f31(I,1),f31(I,3),'-.');plot(f32(I,1),f32(I,3),'-.');
plot(f33(I,1),f33(I,3),'-.');plot(f34(I,1),f34(I,3),'-.');
plot(f35(I,1),f35(I,3),'-.');
axis([0 1 -4 4]);box on
text(0.20,2.4,'5','FontSize',FS);text(0.25,2.2,'4','FontSize',FS)
text(0.33,1.9,'3','FontSize',FS);text(0.43,1.4,'2','FontSize',FS)
text(0.48,0.8,'m=1','FontSize',FS)
xlabel('\itb/a','FontSize',FS);ylabel('\sigma','FontSize',FS)

title('Frequencies buckling modes at neutral stability','FontSize',FS)
text(0.025,-3.65,'Fig. 3 in van Putten, 2002, ApJ, 575, L71','FontSize',12)  

print -deps2 pC.eps; print -djpeg pC.jpg

%--------------------------------------------------------------
% Eqn.(76), van Putten, 2012, Prog. Theor. Phys.,127, 331 
%--------------------------------------------------------------
figure(4);
for k=1:2
for P=1:5
    clear x y;
    if(P==1)delta=f31;N=750;end;if(P==2)delta=f32;N=750;end
    if(P==3)delta=f33;N=200;end;if(P==4)delta=f34;N=200;end
    if(P==5)delta=f35;N=100;end
    for i=1:N
        x(N+i)  =delta(i,1);y(N+i)  =delta(i,2);
        x(N+1-i)=-x(N+i);   y(N+1-i)=y(N+i);
    end
    p=polyfit(x,y,2);PA(P,:)=p;
    if(k==2)plot(delta(:,1),delta(:,2));hold on;end
    if(k==1)a=0:0.01:1;Y=polyval(p,a);plot(a,Y,'*');hold on;end;
end
end
axis([0 1 1.5 2.5]);
legend(['c_{10}=' num2str(PA(1,3)) ', c_{12}=' num2str(PA(1,1))],['c_{20}=' num2str(PA(2,3)) ', c_{22}=' num2str(PA(2,1))],['c_{30}=' num2str(PA(3,3)) ', c_{32}=' num2str(PA(3,1))],['c_{40}=' num2str(PA(4,3)) ', c_{42}=' num2str(PA(4,1))],['c_{50}=' num2str(PA(5,3)) ', c_{52}=' num2str(PA(5,1))],1);

title('Coefficients (c_{m0},c_{m2}) of quadratic fits','FontSize',FS)
text(0.025,1.55,'Eqn. (76) in van Putten, 2012, Prog. Theor. Phys., 127, 331','FontSize',12)  
 
print -deps2 pD.eps; print -djpeg pD.jpg