Hai ragione, scusami, ero convinto che la risoluzione delle immagini fosse più alta. Ma in ogni caso in formato "testo" è meglio così si può anche copiare. I codici sono i seguenti:
Programma principale:
y = linspace (-1, 1, 201);
phi = linspace (-0.016, 0.016, 33);
w0 = 0.5;
d0_in_bs = 50;
alpha = 0;
disp("start");
for i=1:length(y)
disp(y(i))
for j=1:length(phi)
[R_s, w_s, w0_s, x_s, theta_s, R_l, w_l, w0_l, x_l, theta_l] = Archibeams( y(i), phi(j), w0, d0_in_bs, alpha);
C(i,j) = Archicontrast (x_s, w_s, w0_s, x_l, w_l, w0_l);
end
end
[yy, pphi] = meshgrid(y, phi);
mesh(yy, pphi, C);
Funzione "Archibeams":
function [R_s, w_s, w0_s, x_s, theta_s, R_l, w_l, w0_l, x_l, theta_l] = Archibeams( y, phi, w0, d0_in_bs, alpha)
d0_f_bs = 50;
l = 50;
d0_bs_l = 50;
f = 75;
lambda = 633;
lambda = lambda *(10^(-6)); % Per portare la lunghezza d'onda in mm.
% DETERMINAZIONE DI TUTTI I PARAMETRI GEOMETRICI ---------------------------------------------------
% (anche di quelli necessari per la parte gaussiana)
% BRACCIO CORTO:
x_1 = phi*(d0_in_bs + d0_bs_l) + y;
phi_1 = phi - x_1/f;
x_2 = [1/(1 - phi_1*alpha)]*(x_1 + phi_1*(f - alpha*l));
d_ls = f - (l + x_2)*alpha;
x_3 = x_2 - d_ls*(2*alpha - phi_1);
theta_s = 2*alpha - phi_1 + x_3/f;
x_s = x_3 - (d0_bs_l + d0_f_bs)*phi_1;
% BRACCIO LUNGO:
X_1 = phi*(d0_in_bs + 2*l + d0_bs_l) + y;
PHI_1 = phi - X_1/f;
X_2 = [1/(1 + PHI_1*alpha)]*(X_1 + PHI_1*(f + alpha*l));
D_ls = f + (l - X_2)*alpha;
X_3 = X_2 - D_ls*(2*alpha - PHI_1);
theta_l = 2*alpha - PHI_1 + X_3/f;
x_l = X_3 - (d0_bs_l + 2*l + d0_f_bs)*PHI_1;
% DETERMINAZIONE DEI PARAMETRI GAUSSIANI -----------------------------------------------------------
w0_squared = w0*w0;
% BRACCIO CORTO:
z_1 = [(y - x_1)^2 + (d0_in_bs + d0_bs_l)^2]^(1/2);
[R1, w1_squared] = GaussianPropagator_inSpace(w0_squared, z_1, lambda);
[R2, w2_squared, w2_0_squared, z_2] = GaussianPropagator_inLens(R1, w1_squared, lambda, f);
d_lsl = [(x_2 - x_1)^2 + (d_ls)^2]^(1/2) + [(x_3 - x_2)^2 + (d_ls)^2]^(1/2);
[R3, w3_squared] = GaussianPropagator_inSpace(w2_0_squared, z_2 + d_lsl, lambda);
[R4, w4_squared, w4_0_squared, z_4] = GaussianPropagator_inLens(R3, w3_squared, lambda, f);
d_lf = [(x_s - x_3)^2 + (d0_bs_l + d0_f_bs)^2]^(1/2);
[R_s, w_s_squared] = GaussianPropagator_inSpace(w4_0_squared, z_4 + d_lf, lambda);
w0_s = (w4_0_squared)^(1/2);
w_s = (w_s_squared)^(1/2);
% BRACCIO LUNGO:
Z_1 = [(y - X_1)^2 + (d0_in_bs + 2*l + d0_bs_l)^2]^(1/2);
[R1_l, W1_squared] = GaussianPropagator_inSpace(w0_squared, Z_1, lambda);
[R2_l, W2_squared, W2_0_squared, Z_2] = GaussianPropagator_inLens(R1_l, W1_squared, lambda, f);
D_lsl = [(X_2 - X_1)^2 + (D_ls)^2]^(1/2) + [(X_3 - X_2)^2 + (D_ls)^2]^(1/2);
[R3_l, W3_squared] = GaussianPropagator_inSpace(W2_0_squared, Z_2 + D_lsl, lambda);
[R4_l, W4_squared, W4_0_squared, Z_4] = GaussianPropagator_inLens(R3_l, W3_squared, lambda, f);
D_lf = [(x_l - X_3)^2 + (d0_bs_l + 2*l + d0_f_bs)^2]^(1/2);
[R_l, w_l_squared] = GaussianPropagator_inSpace(W4_0_squared, Z_4 + D_lf, lambda);
w0_l = (W4_0_squared)^(1/2);
w_l = (w_l_squared)^(1/2);
Funzione "Archicontrast":
function C = Archicontrast (x_s, w_s, w0_s, x_l, w_l, w0_l)
x_s = x_s*1000;
x_l = x_l*1000;
w_s = w_s*1000;
w0_s = w0_s*1000;
w_l = w_l*1000;
w0_l = w0_l*1000;
x = linspace (-1501, 1500, 3002);
y = linspace (0, 1501, 1502);
[xx_s, yy] = meshgrid((x - x_s).^2 , y.^2);
I1 = sum(sum([(w0_s/w_s)^2]*e.^[-2*(xx_s + yy)/w_s]));
[xx_l, yy] = meshgrid((x - x_l).^2 , y.^2);
I2 = sum(sum([(w0_l/w_l)^2]*e.^[-2*(xx_l + yy)/w_l]));
I3 = sum(sum(2*[(w0_s*w0_l)/(w_s*w_l)]*e.^[-(xx_s + yy)/w_s - (xx_l + yy)/w_l]));
C = I3/(I1 + I2);
Funzione "GaussianPropagator_inSpace":
function [R, w_squared] = GaussianPropagator_inSpace(w0_squared, z, lambda)
lambda = lambda*(10^(-6));
R = z*[1 + ((pi*w0_squared)/(lambda*z))^2];
w_squared = w0_squared*[1 + ((lambda*z)/(pi*w0_squared))^2];
Funzione "GaussianPropagator_inLens":
function [R2, w2_squared, w2_0_squared, z] = GaussianPropagator_inLens(R1, w1_squared, lambda, f)
lambda = lambda*(10^(-6));
R2 = 1/(1/R1 - 1/f);
w2_squared = w1_squared;
w2_0_squared = w2_squared/[1 + ((pi*w2_squared)/(lambda*R2))^2];
z = R2/[1 + ((lambda*R2)/(pi*w2_squared))^2];
Non ho provato col debugger...come si fa per provare?