Homec4science
Diffusion yearII_reports (master) Browse

rSTICAZZI/physnum/rap8/var_n.m
976942d588e1master

Last Change
Raw File

var_n.m
View Options

% Ce script Matlab automatise la production de resultats
% lorsqu'on doit faire une serie de simulations en
% variant un des parametres d'entree.
%
% Il utilise les arguments du programme (voir ConfigFile.h)
% pour remplacer la valeur d'un parametre du fichier d'input
% par la valeur scannee.
%
%% Parametres %%
%%%%%%%%%%%%%%%%
repertoire = './'; % Chemin d'acces au code compile (NB: enlever le ./ sous Windows)
executable = 'exercice8'; % Nom de l'executable (NB: ajouter .exe sous Windows)
input = 'configuration.in'; % Nom du fichier d'entree de base
tfin = 5000;
xl = -200;
xr = 200;
bounds = sprintf('(%.15g, %.15g)', xl, xr);
omega = 0.003;
nsimul = 30;
delt = 64;
sigma = 0.06;
%n = linspace(5, 30, nsimul);
n = linspace(10.9, 11.2, nsimul);
N = 300;
dt = 1;
x0 = -delt;
sampling = 1; %floor(tfin / dt);
hbar = 1;
m = 1;
p0 = hbar * 2 * pi * n / (xr - xl);
A = p0 / (omega * m);
%% Simulations %%
%%%%%%%%%%%%%%%%%
output = cell(1, nsimul);
dP = zeros(1, nsimul);
for i=1:nsimul
output{i} = ['m_outs/nvar_', num2str(n(i))];
cmd = sprintf('%s%s %s N=%d tfin=%.15g dt=%.15g "bounds=%s" omega=%.15g x0=%.15g sigma=%.15g n=%.15g delta=%.15g output=%s sampling=%d outmode=2 >& log', ...
repertoire, executable, input,...
N, tfin, dt, ...
bounds, omega, x0, sigma, n(i), delt, ...
output{i}, sampling);
disp(cmd);
system(cmd);
%psi = load(sprintf('%s_density.out', output{i}));
state = load(sprintf('%s_state.out', output{i}));
%t = state(:,1);
%x = gen_space(xl, xr, N);
ind = optimal_index(state(:,5));
dP(i) = abs(state(ind, 2) -0.5);
end
%% States %%
V0 = potent(m, omega, delt);
figure
plot(n, dP, 'k+')
grid on
xlabel('n')
ylabel('dP')
%% Print values %%
outfile = 'outs/tunnel_var_n.dat';
% print data on file
out = fopen(outfile, 'w');
fprintf(out, 'n dP\n');
for i = 1:nsimul
fprintf(out, '%.15g %.15g\n', n(i), dP(i));
end
fclose(out);
%% Stima %%
estimation = 11.06;
% separe
j = 1;
while(n(j) < estimation)
j = j + 1;
end
nm = n(1:j);
nM = n(j+1:end);
dPm = dP(1:j);
dPM = dP(j+1:end);
fitm = polyfit(nm, dPm, 1);
fitM = polyfit(nM, dPM, 1);
nlimm = - fitm(2) / fitm(1);
nlimM = - fitM(2) / fitM(1);
%% Utils %%
function ind = optimal_index(xmid)
% find max x in time
[~, inds] = findpeaks(smooth(xmid));
ind = inds(1);
end
function V = potent(m, omega, delt)
V = 0.5 * m * omega^2 * delt.^2;
end
function x = gen_space(xl, xr, N)
x = zeros(N+1, 1);
dx = (xr - xl) / N;
for k = 0:N
x(k+1) = k * dx + xl;
end
end
c4science · Help