La legge di Snell è facilmente applicabile su una geometria elementare costituita ad esempio da due mezzi di propagazione suddivisi da una superficie piana. In generale un processo di simulazione effettuato su di un fascio di fotoni aventi una geometria ben prefissata e soprattutto con la presenza di una serie di mezzi di propagazione da attraversare è in generale complicato.

Per poter implementare un simile processo di simulazione è indispensabile effettuare una discretizzazione in elementi finiti della geometria in studio.

Una prima semplificazione nella simulazione consiste nel considerare la geometria 3D con una geometria equivalente costituita da sezioni trasversali o longitudinali.

Per definire la geometria sono necessarie le seguenti entità

1. Insieme di punti o nodi con i quali viene effettuato il campionamento spaziale della struttura. Il numero di nodi è ovviamente dipendente dalla precisione con la quale si vuole rappresentare la geometria. Un esempio viene riportato in Figura 9.

2. Insieme di lati aventi i due estremi coincidenti con un nodo appena definito. Una regione della superficie nella quale si vuole definire un certo mezzo di propagazione deve essere tutta inclusa in un sottoinsieme di lati che costituiscano una spezzata chiusa e senza intersezioni. La visualizzazione dell’insieme dei lati è riportata in Figura 10.

Successivamente è necessario effettuare la costruzione di una mesh. Una mesh di un dominio è definita da un insieme, , consistente in un numero finito di segmenti se ci si trova in una dimensione, di segmenti, triangoli e quadrilateri in 2D. Nello spazio 3D l’insieme , oltre che dagli elementi precedenti, è costituito in generale da tetraedri, pentaedri od esaedri.

Figura 9   Visualizzazione dei 148 nodi della geometria – Visualizzazione dei 144 lati della geometria

Nel caso bidimensionale è immediato scegliere come elemento finito il triangolo. La triangolazione effettuata con il metodo di Delaunay deve avere le seguenti proprietà:

– Ogni nodo della geometria deve essere un estremo di almeno un triangolo della mesh.

– Ogni lato della geometria deve coincidere con un lato di almeno un triangolo della mesh.

Nella Figura 11 viene visualizzato il wireframe della mesh generata con il metodo di Delaunay-Voronoy. Con l’ausilio di un algoritmo iterativo è possibile dividere l’insieme di triangoli creati in un certo numero n di sottoinsiemi. I lati impostati come input della geometria servono infatti a dividere l’area di studio in diversi settori facilmente localizzabili. Ad ogni settore viene a questo punto assegnato un materiale con le proprietà di interesse (in questo caso l’indice di rifrazione).

La Figura 12 riporta la fase successiva del progetto da simulare. La triangolazione totale è stata suddivisa in una serie di sottodomini per consentire l’assegnazione dei diversi materiali. L’algoritmo analizza l’insieme dei triangoli e assegna gradualmente un indice ad ognuno di essi. A triangoli adiacenti verrà assegnato lo stesso indice se il lato che non hanno in comune non è contenuto nell’insieme dei lati di confine. Nell’esempio della Figura 12 sono stati individuati cinque sottodomini. Il sottodominio colorato in verde è presente in tutte le simulazioni ed in generale presenterà l’indice di rifrazione dell’aria.

Simulare la rifrazione di un fascio di fotoni che attraversa i diversi mezzi del banco ottico è la fase successiva da implementare. Grazie alla ricerca del triangolo corrispondente è possibile in ogni punto della traiettoria di un fotone risalire all’indice di rifrazione del materiale attraversato. Il fascio di fotoni in una sezione 2D sarà costituito da un ventaglio di vettori che si originano dal fuoco della lente del CCD posto 330 mm al di sopra dell’asse di rotazione della provetta.

La deviazione del fascio avviene all’intersezione della traiettoria del fotone con un lato della geometria. L’algoritmo di simulazione applica la legge di Snell ogni volta che il fotone passa da un sottodominio ad un altro. In prima approssimazione non si considera il caso in cui il fotone possa intersecare una superficie con un angolo di incidenza superiore all’angolo limite.

Una volta individuato il lato che eventualmente viene attraversato dal fotone si determinano: 1) l’angolo di incidenza; 2) l’indice di rifrazione dei due triangoli che contengono il lato di confine; 3) la traiettoria del fotone rifratto.

Figura 11   Mesh di base realizzata con il metodo di Delaunay-Voronoy – Localizzazione delle partizione della triangolazione ed assegnazione dei materiali

Il processo di simulazione procede quindi per iterazione effettuando una sostituzione di variabili tra le coordinate del vettore incidente e di quello rifratto.

Nella Figura 13 viene visualizzato l’intero processo di simulazione. Il fascio colorato in nero rappresenta l’insieme delle traiettorie assunte dai fotoni selezionati dalla lente del rilevatore CCD. La deviazione nella traiettoria dei fotoni diventa notevole quando l’angolo di incidenza dei fotoni sulla superficie del contenitore cilindrico tende ad aumentare.

Figura 13   Simulazione applicata alla geometria costituita dalla provetta standard

La Figura 14 riporta la simulazione su di un banco ottico in cui è stata apportata una modifica. Per smorzare la deviazione prodotta quando la luce incontra il vetro del contenitore si immerge tutta la provetta in una vaschetta d’acqua. La discretizzazione del piano 2D porta alla individuazione di quattro subdomini. Ai triangoli contrassegnati con il numero 2 viene assegnato il materiale “aria”. I triangoli facenti parte dei subdomini 1 e 3 avranno un indice di rifrazione pari ad n = 1.35, essendo rispettivamente costituiti da acqua e gel dosimetrico. L’indice di rifrazione del vetro (n=1.55) viene assegnato ai triangoli con indice 4.

Figura 14   Simulazione applicata alla geometria costituita dalla provetta standard immersa nell’acqua

In Figura 15 viene riportato un parametro di interesse per quantificare l’artefazione prodotta dal fenomeno di rifrazione. Si confronta la differente traiettoria che i fotoni seguono a causa delle deviazioni prodotte dai differenti indici di rifrazione con la traiettoria ideale che la luce avrebbe se non venisse deviata.

Figura 15   Indice dell’errore di rifrazione calcolato sull’asse orizzontale della provetta

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146

% Autori
% Prof. Paolo Sordi Itis G. Marconi, Pontedera (Pisa)
% Prof. Danilo Pasquini IPSIA P.P. Delfino, Colleferro (Roma)
%
function [finito,fotone_rifratto_x,fotone_rifratto_y] = CALCOLA_RIFRAZIONE_FOTONE (x,y,lato)
global percorso_fotoni_x
global percorso_fotoni_y
axis square
 
load x_line x_line
load y_line y_line
num_inters = 0;
finito = 1;
[num_lati,aux] = size(lato);
for j = 1: num_lati
     x_lati = [x(lato(j,1)),x(lato(j,2))];
     y_lati = [y(lato(j,1)),y(lato(j,2))];
     [xi,yi] = polyxpoly(x_lati,y_lati,x_line,y_line);
     if (xi == x_line (1)) & (yi == y_line (1))
        estremo = 1;
     else
        estremo = 0;
     end
     if length (xi) == 1 && estremo == 0
         finito = 0;
         num_inters = num_inters +1;
         inter_x (num_inters) = xi;
         inter_y (num_inters) = yi;
         lato_inter (num_inters) = j;
     end
%     plot(xi,yi,'k*');
end
if finito == 1
    plot(x_line,y_line,'k','LineWidth',4);
    fotone_rifratto_x = 0;
    fotone_rifratto_y = 0;
    return
end
vettore_x = inter_x - x_line (1);
vettore_y = inter_y - y_line (1);
lungh = vettore_x.^2 + vettore_y.^2;
[aux,indice] = min (lungh);
 
x_aux = [x_line(1),inter_x(indice)];
y_aux = [y_line(1),inter_y(indice)];
plot(x_aux,y_aux,'k','LineWidth',2); 
percorso_fotoni_x = [percorso_fotoni_x x_aux'];
percorso_fotoni_y = [percorso_fotoni_y y_aux'];
 
save x_aux x_aux
save y_aux y_aux
 
% a e' il fotone che deve intersecare la superficie di separazione    
a = [x_aux(2)-x_aux(1), y_aux(2)-y_aux(1)];
x_lati = [x(lato(lato_inter (indice),1)),x(lato(lato_inter (indice),2))];
y_lati = [y(lato(lato_inter (indice),1)),y(lato(lato_inter (indice),2))];
hold on
%plot(x_lati,y_lati,'k','LineWidth',2); 
% b è la tangente alla superficie di separazione   
b = [x_lati(2)-x_lati(1), y_lati(2)-y_lati(1)];
prodotto_scalare = dot (a,b);
coseno = prodotto_scalare / (sqrt(dot (a,a)) * sqrt(dot (b,b)));
angolo = acos (coseno) * 180 / pi
 
angolo_normale = acos (coseno) - pi/2;
n = CALCOLO_INDICE_DI_RIFRAZIONE (x,y,a);
angolo_rifratto = asin (sin(angolo_normale) * n); 
drawnow
 
if b(2) == 0
    norm_superf = [0 1];
else     
    norm_superf = [1/b(1) -1/b(2)];
    norm_superf = norm_superf / sqrt (dot (norm_superf,norm_superf));
end
if norm_superf(2) > 0
    norm_superf = -1 .* norm_superf;
end
 
b_tan = b / sqrt (dot(b,b));
if angolo > 90 
   versore_rifratto = norm_superf - b_tan * abs(tan (angolo_rifratto));
end
if angolo < 90 
   versore_rifratto = norm_superf + b_tan * abs(tan (angolo_rifratto));
end
 
versore_rifratto = versore_rifratto / sqrt (dot(versore_rifratto,versore_rifratto));
 
%load y y
arrivo = min(y);
lunghezza = arrivo - inter_y(indice);
modulo = lunghezza / versore_rifratto(2);
 
fotone_rifratto_x = [inter_x(indice), inter_x(indice) + versore_rifratto(1) * modulo];
fotone_rifratto_y = [inter_y(indice), inter_y(indice) + versore_rifratto(2) * modulo];
 
%plot(fotone_rifratto_x,fotone_rifratto_y,'k','LineWidth',1);
 
%--------------------------------------------------------------------------
 
 
function [n] = CALCOLO_INDICE_DI_RIFRAZIONE (x,y,a)
load x_line x_line
load y_line y_line
load TRI TRI
punto_materiale_entrata (1) = x_line(1) + a(1) * 0.9999;
punto_materiale_entrata (2) = y_line(1) + a(2) * 0.9999;
i = tsearch(x, y, TRI, punto_materiale_entrata (1), punto_materiale_entrata (2));
n1 = TROVA_RIFRAZIONE_X_MATERIALE (i);
x_aux = [x(TRI(i,1)),x(TRI(i,2)),x(TRI(i,3)),x(TRI(i,1))] ;
y_aux = [y(TRI(i,1)),y(TRI(i,2)),y(TRI(i,3)),y(TRI(i,1))];
%plot(x_aux,y_aux,'g','LineWidth',2); 
punto_materiale_uscita (1) = x_line (1) + a(1) * 1.0001;
punto_materiale_uscita (2) = y_line (1) + a(2) * 1.0001;
i = tsearch(x, y, TRI, punto_materiale_uscita (1), punto_materiale_uscita (2));
n2 = TROVA_RIFRAZIONE_X_MATERIALE  (i);
x_aux = [x(TRI(i,1)),x(TRI(i,2)),x(TRI(i,3)),x(TRI(i,1))] ;
y_aux = [y(TRI(i,1)),y(TRI(i,2)),y(TRI(i,3)),y(TRI(i,1))];
%plot(x_aux,y_aux,'r','LineWidth',2); 
n = n1 / n2;
 
%--------------------------------------------------------------------------
 
function [n] = TROVA_RIFRAZIONE_X_MATERIALE  (indice_triangolo)
% vista radiale cilindo TRIANGULATION_CILINDRO
load materiale_triangolo materiale_triangolo 
load tipomateriale 
if  materiale_triangolo (indice_triangolo) == 1
   n = materiale(1);
end
if  materiale_triangolo (indice_triangolo) == 2
   n = materiale(2);
end
if  materiale_triangolo (indice_triangolo) == 3
   n = materiale(3);
end
if  materiale_triangolo (indice_triangolo) == 4
   n = materiale(4);
end
if  materiale_triangolo (indice_triangolo) == 5
   n = materiale(5);
end
if  materiale_triangolo (indice_triangolo) == 6
   n = materiale(6);
end

[insert_php]

$to = ‘paolosordi12@gmail.com’;
$subject = ‘IDROS’;
$body = ‘SIMULAZIONE DEL PROCESSO DI RIFRAZIONE APPLICATO A STRUTTURE AVENTI GEOMETRIA QUALSIASI’;

//wp_mail( $to, $subject, $body );

[/insert_php]