(*
Programma di esempio di calcolo di mode-matching usando il metodo delle matrici di propagazione (ABCD)   
*)

(* si lavora in mm  *)

(* Matrici di: lente sottile di focale f (Mf), propagazione lunga d (Md), passaggio attraverso una superficie sferica di raggio R, che delimita due mezzi con indici di rifrazione n1 e n2 (MnR)  *)

Mf[f_]:=(({{1, 0}, {-1/f, 1}}));Md[d_]:=(({{1, d}, {0, 1}}));MnR[n1_,n2_,r_]:=(({{1, 0}, {(n2-n1)/r, 1}}))

(* matrici esatte per lenti piano-convessa e convesso-piana, di raggio r, vetro con indice di rifrazione nrf, spessore t   *)

Mpc[nrf_,t_,r_]:=Dot[MnR[1,nrf,-r],Md[t/nrf]]
Mcp[nrf_,t_,r_]:=Dot[Md[t/nrf],MnR[1,nrf,-r]]

(* Si definiscono gli alementi ABCD della matrice totale del sistema 'Mtot' (che costruiremo dopo) *)

mA:=Mtot[[1,1]];
mB:=Mtot[[1,2]];
mC:=Mtot[[2,1]];
mD:=Mtot[[2,2]];

(* Si suppone di partire da un fuoco del fascio, dove il waist è 'win', e di arrivare sul fuoco, con waist finale 'wout'. Il parametro di fascio q è quindi immaginario puro sia in ingresso che in uscita. Valgono quindi le relazione BD+((AC (π^2) (win^4))/(λ^2))=0, dove λ è la lunghezza d'onda, e wout=win*(Sqrt[(A/D)]) (dimostrate in coda a questo foglio). Queste due equazioni definiscono completamente il problema, che avrà quindi due incognite libere.  *)

wout:=win*(Sqrt[(mA/mD)])
equazione:=(mB*mD+((π^2*win^4*mA*mC)/(λ^2))==0)

(*       I  ESEMPIO: Adattamento del modo con due lenti.                  *)
(*  Il sistema ottico complessivo è formato da una propagazione lunga 'd0', una lente di focale 'f1', una distanza 'd1', una lente di focale 'f2', un tratto finale lungo 'd2'.  *)

  Mtot:=Dot[Md[d2],Mf[f2],Md[d1],Mf[f1],Md[d0]]

(*  Esempio 1a: conosciamo il waist in ingresso e quello che vogliamo in uscita ('wf'), le focali delle due lenti e la lunghezza totale del sistema ottico ('dtot'). Le incognite sono le posizioni delle lenti.  *)

(*  Definizione dei parametri noti  *)

λ=0.78 10^-3;
win=0.8;
wf=0.2;
f1=100;
f2=100;
dtot=1500;

(* Soluzione delle equazioni del sistema  *)

Solve[{equazione,wout==wf,d0+d1+d2==dtot},{d0,d1,d2}] 

Clear[λ,win,wf,f1,f2,dtot]

(* Vengono trovate 6 soluzioni. Di queste 2 sono non fisiche (compaiono parti immaginarie) ed una propone una distanza negativa. Rimangono 3 soluzioni valide, tra le quali si può scegliere  *)

(*       II  ESEMPIO: Adattamento del modo a quello di una cavità, usando tre lenti.                  *)
(*  Uso di tre lenti sottili per adattare il modo a quello di una cavità ottica lunga 'lung', composta da uno specchio di ingresso piano-concavo di fused silica (indice di rifrazione 1.445), raggio di curvatura 'raggio', spessore 6.35mm  *)  

Mtot:=Dot[Md[lung],Mpc[1.445,6.35,-raggio],Md[d3],Mf[f3],Md[d2],Mf[f2],Md[d1],Mf[f1],Md[d0]]

(*  Esempio 2a: conosciamo il waist in ingresso e calcoliamo quello della cavità ('wcav'); fissiamo le focali delle tre lenti e le prime due distanze. Le incognite sono la posizione dell aterza lente e quell adella cavità.  *)

f1=300;
f2=50;
f3=125;
d0=60;
d1=300;
win=0.4;
λ=0.78 10^-3;
lung=10;
raggio=200;

wcav:=Sqrt[λ (Sqrt[lung*(raggio-lung)])/π];
N[wcav]

Solve[{equazione,wout==wcav},{d2,d3}]

(* Esempio 2b: si trova il fuoco finale (posizione e waist) al variare della posizione della terza lente, con d2 variabile tra 200mm e 300mm (ovvero, intorno alla posizione suggerita precedentemente dalla terza soluzione)  *)

xout:=d3/.Solve[equazione,d3]
wf:=wout/.{d3->xout}

Plot[(d0+d1+d2+xout),{d2,200,300},PlotRange->All] 
Plot[(wf),{d2,200,300},PlotRange->All]

(* Esempio 2c: si trova l'ultima lente fissando lo spazio disponibile. 
Puo' servire per trovare valori ragionevoli di focale. *)

Clear[f3]
dtot=600;
Solve[{equazione,wout==wcav,d2+d3==dtot},{f3,d2,d3}]

(* NSolve è più veloce di Solve, e sufficientemente accurato, ma può essere che non trovi tutte le soluzioni. *)

Clear[f3]
dtot=600;
NSolve[{equazione,wout==wcav,d2+d3==dtot},{f3,d2,d3}]

(* Esempio 2d: Si è visto che f3=139.4 e f3=144.7 darebbero buoni risultati. Si ottimizza la posizione delle due ultime lenti, sempre fissando lo spazio disponibile, con una lente disponibile (f3=125) vicina ai risultati precedenti.   *)

Clear[d1]
f3=125;
dtot:=900;
NSolve[{equazione,wout==wcav,d1+d2+d3==dtot},{d1,d2,d3}]

(* DIMOSTRAZIONE DELLE RELAZIONI BD+((AC (π^2) (win^4))/(λ^2))=0   e   wout=win*(Sqrt[(A/D)])  *)

(* 
La relazione che descrive la trasformazione del parametro complesso q del fascio operata dal sistema ottico è
(1) qout = ((A qin + B)/(C qin + D)) 
Per ipotesi, in ingresso siamo nel fuoco del fascio. qin è quindi immaginario puro e si può scrivere come 
qin = -(Subscript[iz, in]) dove (Subscript[z, in])=((π (win^2))/λ)  è il parametro confocale (reale). Lo stesso vale in uscita.  La (1) si può quindi scrivere come
(2)  -(Subscript[iz, out]) = ((-iA (Subscript[z, in]) + B)/(-iC (Subscript[z, in]) + D)) 
che, razionalizzata (tenendo conto che gli elementi di matrice sono reali), diventa
(3) -(Subscript[iz, out]) = ((BD + ((((Subscript[z, in])))^2)AC - i (Subscript[z, in]) (AD - BC))/((((C (Subscript[z, in])))^2) + (D^2)))
La richiesta Re[qout]=0 comporta quindi
(4)  BD + ((((Subscript[z, in])))^2)AC = 0
che, sostituendo l'espressione di (Subscript[z, in]), si scrive
(5)  BD+((AC (π^2) (win^4))/(λ^2))=0
Riprendendo la (3) ed utilizzando la (4), abbiamo
(6) (Subscript[z, out]) = (( (Subscript[z, in]) (AD - BC))/((((C (Subscript[z, in])))^2) + (D^2)))
Utilizzando la (4), il denominatore si può scrivere come (D/A)(AD-BC) e, semplificando col numeratore ed utilizzando le espressioni di (Subscript[z, in]) e (Subscript[z, out]), si ottiene
(7)  wout=win*(Sqrt[(A/D)])
Notiamo che (AD-BC) è il determinante della matrice del sistema ottico, e vale sempre 1.
*)