%arrival time of a radio pulse as a function of frequency, nu
%for cold plasma, arrival time is propto nu^p where p=-2 when
%nu>>nu_p, the plasma frequency
%This program determines p as a funciton of nu/nu_p
%
    measured_p=2.006;	%reported arrival time index

    req_index=measured_p/2;
    
	%vg is the group velocity where x=(omega_p/omega)^2
    vg=@(x)((1-x).^(1/2));
	%ivg is 1/vg and is a proxy for arrival time
    ivg = @(x)((1-x).^(-1/2)-1.0); 
    
    epsilon_min=1E-14;
    epsilon_max=1-epsilon_min;
    nsteps=1E5;
    x=linspace(epsilon_min,epsilon_max,nsteps);
    yvg=vg(x);
    yivg=ivg(x);
%%
figure(1)
    subplot(3,1,1)
        semilogy(x,yivg,'-k','LineWidth',1)
        xlim([1e-5 Inf])
        ylim([1e-5 1])
         ylabel('t_a')
         xlabel('\epsilon')
         set(gca,'FontSize',16,'LineWidth',2)

%%
    x1=x; x1(end)=[];
    x2=x; x2(1)=[];
    y1=yivg; y1(end)=[];
    y2=yivg; y2(1)=[];
    slope=log(y1./y2)./log(x1./x2);

    subplot(3,1,2)
    plot(x1,slope,'-k','LineWidth',1)
    ylim([1 1.5])
    xlabel('\epsilon')
    ylabel('{\it q}')
    set(gca,'FontSize',16,'LineWidth',2)
    
%%
    

    subplot(3,1,3)
    plot(x1,slope,'-k','LineWidth',1)
    ylim([1 1.007])
    hold on
    k=find(slope>req_index,1,'first');
    fprintf('slope>%d epsilon=%0.5g\n',req_index,x1(k));
    plot(x1(k),slope(k),'+r','MarkerSize',13)
    xlabel('\epsilon')
    ylabel('{\it q}')
    set(gca,'FontSize',16,'LineWidth',2)
    hold off
    print('ta-nu-index','-dpng')

    
