growthrate.m
0b252c046024master

growthrate.m
View Options

			resfile='/misc/nonlinear/results_testl6_high.h5';

			M=nonlinhdf5(resfile);
			t=M.t2d;

			%% Analytical growth rate
			l=6; % mode number
			gamma_ana=analytic_growth(l,M.ra,M.rb,M.phia,M.phib,M.B0,M.n0,M.rm,M.rp)

			if 1
			%% Growth rate determined with the potential perturbation
			phi=M.phil(:,:,:); % modes
			maxphil=zeros(size(t));
			for n=1:length(t)
			maxphil(n)=max(phi(:,l+1,n));
			end

			% Fit to determine the growth rate
			teval=500:2000;
			p1=polyfit(t(teval),log(maxphil(teval)),1);
			gamma_num1=p1(1)

			figure
			plot(t,log(maxphil))
			hold on
			plot(t,polyval(p1,t),'--')
			xlabel('t [s]')
			ylabel('log \|\phi_l\| [a.u.]')
			title(sprintf('\\gamma_{ana} = %1.3g s^{-1}, \\gamma_{num} = %1.3g s^{-1}',gamma_ana,gamma_num1))
			end

			if 1
			%% Growth rate determined with the electrostatic energy
			Wthet=M.Wthet;

			% Fit to determine the growth rate
			teval=2000:4000;
			p2=polyfit(t(teval),log(Wthet(teval)),1);
			gamma_num2=p2(1)/2

			figure
			plot(t,log(Wthet))
			hold on
			plot(t,polyval(p2,t),'--')
			xlabel('t [s]')
			ylabel('log W_{\theta} [a.u.]')
			title(sprintf('\\gamma_{ana} = %1.3g s^{-1}, \\gamma_{num} = %1.3g s^{-1}',gamma_ana,gamma_num2))
			end