Università degli Studi di Roma “Tor Vergata”
FACOLTÀ DI SCIENZE MATEMATICHE FISICHE E NATURALI
Corso di Laurea Specialistica in Scienze dell'Universo
TESI DI LAUREA
Ottimizzazione della ricostruzione del fronte d'onda
per telescopi solari di grande apertura
Candidato:
Marco Stangalini
Matricola 0103666
Relatore:
Prof. Francesco Berrilli
Relatore esterno:
Dott. Dario Del Moro
ANNO ACCADEMICO 2007-2008
La turbolenza atmosferica produce rapide fluttuazioni dell’indice di rifrazione
dell’aria che contribuiscono in gran misura alla degradazione della risoluzione
spaziale delle immagini riprese dai telescopi di Terra. L’Ottica Adattiva è
una tecnologia in grado di porre rimedio a questo inconveniente effettuando
una compensazione in tempo reale degli effetti della turbolenza atmosferica.
Negli ultimi anni, i progressi dell’Astrofisica hanno reso necessario lo sviluppo di questa tecnologia al fine di ottenere immagini degli oggetti celesti con
risoluzione spaziale confrontabile, o superiore, a quella raggiunta dallo spazio.
Un capitolo importante dell’Ottica Adattiva è quello associato alle osservazioni solari. In questo caso l’effetto del seeing indotto dall’irraggiamento
solare diurno e la presenza di immagini non puntiformi, richiedono lo sviluppo di particolari accorgimenti. Le regioni esterne del Sole sono altamente
strutturate e dinamiche su una vasta gamma di scale spaziali. Capire la
fisica delle strutture solari osservate su piccola scala rappresenta un requisito fondamentale alla comprensione di moltissimi aspetti ancora poco chiari:
l’immagazzinamento dell’energia magnetica su piccola scala e il suo ruolo nei
fenomeni coronali, la generazione e la dissipazione del campo magnetico e il
riscaldamento coronale. In questo quadro, l’ottica adattiva diviene una tecnologia molto importante per i moderni telescopi solari di grande apertura,
e i moderni sviluppi di questa disciplina sembrano essere molto promettenti.
In questo lavoro di tesi, si è lavorato principalmente sull’ottimizzazione del
processo di ricostruzione del fronte d’onda. Partendo dai dati acquisiti presso
il telescopio solare Franco-Italiano THEMIS, si è cercato il miglior compromesso tra tempo di calcolo e accuratezza della descrizione delle aberrazioni,
cercando di fornire risultati in grado di indirizzare gli sforzi verso l’uso di tecniche di ricostruzione basate su metodi PCA. In un secondo momento è stato
implementato un codice di simulazione dell’effetto delle turbolenze, adatto allo studio del comportamento di telescopi solari di grande diametro. Attraverso queste simulazioni è stato possibile valutare la qualità della correzione
adattiva facendo uso di tecniche statistiche per la riduzione dell’errore associato al ritardo temporale dei sistemi di compensazione. Tali tecniche, basate
su modelli autoregressivi di sequenze temporali, sembrano essere un prezioso
ausilio verso la correzione ad alta frequenza delle turbolenze atmosferiche.
Indice
1 Principi di ottica adattiva solare
1.1 Fluttuazioni dell’indice di rifrazione dell’aria: le aberrazioni
dinamiche . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Distorsioni di fase nel fronte d’onda . . . . . . . . . . . . . .
1.3 Ottica Adattiva . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 PSF di un sistema ottico . . . . . . . . . . . . . . . . . . . .
1.5 Il Seeing diurno . . . . . . . . . . . . . . . . . . . . . . . . .
1.6 I sensori di fronte d’onda . . . . . . . . . . . . . . . . . . . .
1.7 L’algoritmo SAD . . . . . . . . . . . . . . . . . . . . . . . .
1.8 Gli specchi deformabili . . . . . . . . . . . . . . . . . . . . .
1.9 Schema di controllo e chiusura del loop di correzione . . . .
3
.
.
.
.
.
.
.
.
.
4
5
8
9
9
11
14
14
14
2 Caratterizzazione modale del fronte d’onda
2.1 Espansione del fronte d’onda in polinomi di Zernike . . . . . .
2.2 Ricostruzione della fase . . . . . . . . . . . . . . . . . . . . . .
2.3 PCA: i polinomi di Karhunen-Loeve . . . . . . . . . . . . . . .
2.4 Ricostruzione modale per il telescopio solare THEMIS . . . . .
2.4.1 Contributo al Tip-Tilt dei modi di alto ordine . . . . .
2.4.2 Algoritmo IDL per la ricostruzione modale di KarhunenLoeve . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5 Fitting error . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5.1 L’errore di campionamento . . . . . . . . . . . . . . . .
2.6 Conclusioni . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
21
24
25
27
28
3 Simulazioni di turbolenza multi-layer
3.1 Turbolenza . . . . . . . . . . . . . . . . . . .
3.2 Simulazioni basate su FFT . . . . . . . . . . .
3.3 Implementazione multi-layer . . . . . . . . . .
3.3.1 L’algoritmo di simulazione multi-layer
3.4 Conclusioni . . . . . . . . . . . . . . . . . . .
37
37
41
44
45
47
3
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
28
30
32
33
4 Ottimizzazione della sequenza di correzione attraverso previsioni di turbolenza basate su modelli ARMA
4.1 I processi stocastici . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 Stazionarietà . . . . . . . . . . . . . . . . . . . . . . .
4.1.2 Processi ARMA . . . . . . . . . . . . . . . . . . . . . .
4.2 Previsione di un processo stazionario nel dominio del tempo .
4.3 Previsione dei coefficienti di Zernike . . . . . . . . . . . . . . .
4.4 Conclusioni . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5 Studio della stabilità del processo di compensazione attraverso simulazioni di controreazione
5.1 Il sistema AO come sistema controreazionato . . . . . . . . . .
5.2 Simulazione del loop di controreazione . . . . . . . . . . . . .
5.3 Stabilità del sistema SLS in funzione del numero di polinomi
utilizzati per la ricostruzione . . . . . . . . . . . . . . . . . . .
5.4 Conclusioni . . . . . . . . . . . . . . . . . . . . . . . . . . . .
49
50
50
52
55
58
64
65
66
67
69
73
Conclusioni generali
75
Appendice sorgenti IDL
81
Pubblicazioni
121
Lista acronimi utilizzati
• ACF Auto Correlation Function
• ACVF Auto Covariance Function
• AO Adaptive Optics
• ARMA Auto Regressive Moving Average processes
• DM Deformable Mirror
• EST European Solar Telescope
• FWHM Full Width Half Maximum
• PCA Principal Component Analisys
• PSF Point Spread Function
• WFS Wavefront sensor
• RTC Real Time Computer
• SAD Sum of Absolute Differences
• SLE Servo Lag Error
• SLS Servo Loop System
• SVD Singular Value Decomposition
• TDE Time Delay Error
• THEMIS Télescope Héliographique pour l’Étude du Magnétisme et des
Instabilités Solaires
1
2
Capitolo 1
Principi di ottica adattiva
solare
L’atmosfera solare è altamente strutturata e dinamica. Capire la fisica delle
strutture su piccola scala è in molti casi di cruciale importanza per fornire
una risposta ai grandi problemi, ancora aperti, che la fisica solare ci riserva.
Modelli teorici di tubi di flusso, ad esempio, prevedono strutture magnetiche
su scale spaziali di circa 70 Km sulla fotosfera solare, ben al di sotto dell’attuale potere risolutivo dei telescopi solari. L’osservazione su queste scale
spaziali è perciò di fondamentale importanza, ma richiede risoluzioni angolari
di almeno 0.1 arcosecondi (circa 70 Km sulla superficie solare). Sebbene sia
possibile ottenere immagini al limite diffrattivo con l’uso di tecniche di restoring post-facto, molto spesso per studi spettropolarimetrici e spettroscopici è
necessario avere, allo stesso tempo, alta risoluzione spaziale e spettrale.
Una buona soluzione in grado di fornire una risposta a queste richieste è
rappresentata dall’uso di un sistema di compensazione adattiva. La progettazione e la realizzazione di un sistema di Ottica Adattiva (AO) richiede
lo studio e l’ottimizzazione della ricostruzione delle aberrazioni presenti nel
fronte d’onda incidente.
L’impiego dell’Ottica Adattiva sui telescopi solari, a differenza di quelli notturni, risulta particolarmente complessa a causa della presenza di un target
a basso contrasto ed esteso, quale appunto la granulazione fotosferica solare.
Da questi due soli fattori scaturiscono una vasta gamma di problematiche
che richiedono uno studio approfondito nella realizzazione dei sistemi di AO
solare.
In questo capitolo vedremo una breve introduzione ai principali concetti dell’ottica adattiva, con un particolare accento sulle problematiche connesse
all’uso di questa tecnologia nelle osservazioni solari.
3
1.1
Fluttuazioni dell’indice di rifrazione dell’aria: le aberrazioni dinamiche
Qualsiasi sistema ottico reale possiede delle imperfezioni tali da contribuire
alla non parassialità dei raggi luminosi. Queste imperfezioni contribuiscono
alle aberrazioni intrinseche del sistema ottico. In completa analogia si può
dire che le fluttuazioni dinamiche dell’indice di rifrazione del mezzo in cui si
propaga un onda luminosa determinano piccoli scarti dalla parassialità dei
raggi, con la conseguente comparsa di aberrazioni nelle immagini risultanti.
Le fluttuazioni dell’indice di rifrazione del mezzo possono essere molto rapide
e sono essenzialmente proporzionali alle fluttuazioni puntuali di temperatura. Il loro effetto è perciò altamente dinamico su una vasta gamma di scale
spaziali e temporali.
L’atmosfera terrestre può essere pensata come una sovrapposizione di strati
turbolenti con delle fluttuazioni dell’indice di rifrazione locale, la cui statistica segue quella delle fluttuazioni di temperatura descritte dal modello di
turbolenza di Kolmogorov.
Supponiamo di avere in un certo punto r un indice di rifrazione n(r) e un
indice n(r + ρ) nel punto (r + ρ). La varianza della differenza dei due
valori è data da una funzione di struttura dell’indice di rifrazione, definita
nel modo seguente:
2
Dn (ρ) = h|n(r) − n(r + ρ)|2 i = CN2 ρ 3
(1.1)
dove i simboli h...i rappresentano la media sull’ensemble e CN2 il coefficiente di
struttura, il quale è una misura delle disomogeneità locali dell’indice stesso.
In prima approssimazione, la funzione di struttura dell’indice dipende solo
dal modulo della distanza tra i due punti e non dalla loro posizione, questo
ci suggerisce che il processo casuale sia, almeno localmente, omogeneo. La
relazione (1.1) risulta valida solo per valori di ρ più piccoli della ’outer scale’
di turbolenza, ossia la scala spaziale in cui viene iniettata l’energia e da cui
parte la cascata dissipativa verso scale sempre più piccole. Su scale spaziali
più grandi le fluttuazioni diventano scorrelate e la funzione di struttura non
dipenderà più da ρ tendendo quindi ad infinito e rendendo la trattazione
inadeguata.
La scelta della scala di iniezione di energia è stata ampiamente dibattuta
in letteratura, la gamma di valori proposti va da poche decine di centimetri
a diversi Km, in base al sistema fisico in esame. Nei capitoli successivi,
affrontando il problema della simulazione delle turbolenze atmosferiche, vedremo che l’atmosfera può essere schematizzata come la sovrapposizione di
diversi layer, ognuno con la propria scala caratteristica dettata dalla quo4
ta dello strato turbolento. Sperimentalmente è stato appurato che la legge
di potenza (1.1) risulta accurata su distanze minori di un metro (Roddier, 2004). L’equazione precedente descrive statisticamente la distribuzione
spaziale delle fluttuazioni ad un certo istante di tempo t. Certamente le fluttuazioni avranno una dipendenza temporale oltre che spaziale e per questo
motivo si può, in analogia, definire una funzione di struttura che descriva
l’evoluzione temporale delle fluttuazioni lungo la linea di vista:
DN (τ ) = h|n(r, t) − n(r, t + τ )|2 i
(1.2)
Il tempo di vita di una fluttuazione dell’indice di rifrazione è sperimentalmente di gran lunga più lungo del tempo necessario all’attraversamento della
linea di vista da parte della fluttuazione stessa. Ciò è vero per una turbolenza ’open air ’ e perturbazioni guidate dal vento. Questa approssimazione è
chiamata approssimazione di Taylor o del flusso congelato:
n(r, t + τ ) = n(r − vτ, t)
(1.3)
dove v è la velocità del vento che “trasporta“ la fluttuazione.
Introducendo la (1.3) nella (1.1) si ottiene la seguente espressione per l’andamento della funzione di struttura temporale:
DN (τ ) = h|n(r, t) − n(r − vτ ), t|2 i = CN2 |vτ |2/3
(1.4)
la funzione di struttura temporale si ottiene, perciò, sostituendo |vτ | a ρ
nella funzione di struttura spaziale.
L’indice di rifrazione del mezzo non dipende solo dalla temperatura ma
anche dalla lunghezza d’onda di osservazione, per questo motivo anche le
funzioni di struttura dipenderanno a loro volta da questa quantità.
1.2
Distorsioni di fase nel fronte d’onda
Le onde possono essere convenientemente descritte attraverso una funzione
complessa Ψ:
Ψ = A · exp(iϕ)
(1.5)
dove A e ϕ sono numeri reali che rappresentano rispettivamente l’ampiezza e la fase del campo di fluttuazione. Una superficie con lo stesso Ψ è
chiamata fronte d’onda. Prima di incontrare l’atmosfera il fronte d’onda
proveniente dalle sorgenti astronomiche è piatto. All’interno dell’atmosfera
la velocità della luce varierà come l’inverso dell’indice di rifrazione. I raggi
5
luminosi che si propagano in zone di alto indice di rifrazione subiranno un ritardo di fase maggiore rispetto agli altri. Il fronte d’onda risultante sarà una
superficie non più piatta bensı̀ corrugata. Come accennato in precedenza, i
ritardi di fase non sono acromatici. A maggiori lunghezze d’onda corrispondono, infatti, fluttuazioni più piccole.
La deformazione del fronte d’onda è legata alle fluttuazioni del cammino
ottico attraverso la relazione:
Z
δ = n(z)dz
(1.6)
Il principio della AO è che tali ritardi di fase possono essere compensati
attraverso uno specchio deformabile che assuma la stessa forma del fronte
d’onda, con la metà dell’ampiezza delle fluttuazioni.
I ritardi di fase possono essere messi in relazione alla variazione dell’indice
di rifrazione:
Z
ϕ = k n(z)dz
(1.7)
dove k è il numero d’onda (k = 2π/λ) e n(z) l’indice di rifrazione lungo
la linea di vista descritta dalla variabile z. E’ evidente, a questo punto,
la dipendenza del ritardo di fase dalla lunghezza d’onda della radiazione
incidente. Ciò che più conta non è il valore assoluto della fase in ciascun punto
del fronte d’onda ma la differenza di fase tra due punti separati spazialmente
da un vettore ξ.
Se ϕ è la fase del fronte d’onda, possiamo quindi definire un funzione di
struttura di fase come:
Dϕ (ξ) = h|ϕ(x) − ϕ(x + ξ)|2 i
(1.8)
Usando le relazioni precedenti si può scrivere:
Z
5
2
Cn2 (z)dzξ 3
Dϕ (ξ) = 2.91k
(1.9)
dove l’integrale è calcolato lungo la linea di vista. L’atmosfera è generalmente
considerata stratificata in piani paralleli, ciò significa che CN2 dipende solo
dalla quota h dal suolo:
Z
5
2
−1
Cn2 (h)dhξ 3
Dϕ (ξ) = 2.91k (cosγ)
(1.10)
dove cosγ rappresenta la massa d’aria. Talvolta risulta molto utile riscrivere
questa relazione nel seguente modo:
Dϕ (ξ) = 6.88(
6
ξ 5
)3
r0
(1.11)
dove:
2
r0 = [0.423k (cos(γ))
−1
Z
5
CN2 (h)dh]− 3
(1.12)
è chiamato parametro di Fried ed ha le dimensioni di una lunghezza spaziale.
Questo parametro fornisce la misura delle dimensioni entro le quali una certa
zona del fronte d’onda risulta isoplanatica1 .
E’ utile notare come il parametro di Fried cresca all’aumentare della lunghezza d’onda.
In base alle definizioni di FRIED (1965) e Noll (1976):
5
σ 2 = 1.03(d/r0 ) 3
(1.13)
perciò, se si considera una certa area del fronte d’onda di raggio d pari a r0 , lo
sfasamento rms sulla circonferenza risulta essere pari ad un radiante. Questa
può essere presa in considerazione come definizione stessa del parametro r0 .
L’equazione 1.11 descrive la distribuzione spaziale delle distorsioni del fronte
d’onda. Essa ci consente di determinare i parametri necessari a descrivere
l’ampiezza della correzione da applicare nei sistemi AO, ma non ci dice quanto velocemente ciò deve essere fatto. Per risolvere questo problema si può
definire, in analogia a quanto fatto in precedenza, una funzione di struttura
di fase dipendente dal tempo:
Dϕ (τ ) = h|ϕ(x, t) − ϕ(x, t + τ )|2 i
(1.14)
essa ci fornisce indicazioni sulla variabilità temporale dei ritardi di fase e sulla loro correlazione. A questo punto si può associare al fronte d’onda anche
un tempo di coerenza τ0 , in aggiunta al parametro di Fried che definisce la
coerenza spaziale di un patch isoplanatico.
Se infatti consideriamo l’atmosfera come una stratificazione di diversi
layer di turbolenza che si muvono con velocità diverse, possiamo definire un
velocità media pesata sui diversi layer v0 , la quale ci fornisce:
τ0 = r0 /v0
(1.15)
E’ bene sottolineare come questa approssimazione, e con essa quella di Taylor,
non prendano in esame la presenza di shear e vorticità nei flussi turbolenti.
1
Con isoplanatico si intende quella regione del fronte d’onda i cui estremi sono sfasati
per una quantità minore di un radiante.
7
Figura 1.1: Il fronte d’onda piatto proveniente dalla sorgente nell’incontrare l’atmosfera
terrestre viene corrugato a causa della turbolenza.
1.3
Ottica Adattiva
Come accennato brevemente in precedenza, i ritardi di fase possono essere
compensati mediante l’uso di uno specchio deformabile, una volta descritto
il fronte d’onda incidente. Quando si progetta un sistema AO, è necessario
comprendere quale relazione ci sia tra il fronte d’onda aberrato nell’apertura
del telescopio e la distribuzione delle intensità sul piano focale. Questa relazione è generalmente non lineare. Un incremento lineare nell’ampiezza di
queste perturbazioni produce effetti esponenzialmente crescenti nell’immagine risultante (Roddier, 2004).
Abbiamo già visto che la distorsione di fase rms è per definizione pari ad
un radiante in una circonferenza di raggio r0 . Su un telescopio di diametro
D il numero di queste aree di coerenza è pari a D/r02. Questo è un ordine
di grandezza per il numero di parametri da controllare. Se si utilizza uno
specchio segmentato, ad esempio, si può dimostrare (Roddier, 2004) che il
numero minimo di attuatori in grado di riportare una fluttuazione di fase al
di sotto di un radiante è:
D
N = 0.27( )2
(1.16)
r0
Dato che r0 ∝ λ6/5 , il numero di attuatori decresce come λ−12/5 . Ciò sottolinea ancora una volta che è molto più semplice compensare i ritardi di fase
8
a lunghezze d’onda grandi. Se si considera per r0 un valore tipico di 15cm
e D = 4m, come lo specchio previsto per il telescopio europeo EST, si può
stimare il numero minimo di attuatori necessari alla correzione che risulta
essere pari a 192. Si può anche definire un ritardo temporale τ0 dettato dalla
velocità del vento:
τ0 = 0.314 · (r0 /v)
(1.17)
questo ritardo è chiamato ritardo di Greenwood (Fried, 1990).
1.4
PSF di un sistema ottico
Ogni sistema ottico ha una propria funzione di risposta ad un impulso luminoso. Tale funzione è detta PSF (Point Spread Function) e rappresenta
il pattern di diffrazione prodotto in condizioni ideali dal fronte d’onda nella
propagazione all’interno del sistema ottico. Ciò è dovuto al fatto che qualsiasi sistema, avendo dimensioni finite, riesce a collezionare solo una piccola
porzione del fronte d’onda. L’immagine di una sorgente risulta pertanto la
convoluzione dell’oggetto con la PSF dello strumento.
La PSF, come mostrato in figura 1.2, risente tuttavia della presenza di aberrazioni. In condizioni ideali la PSF è descritta da una funzione di Airy (per
pupille circolari) la cui ampiezza a metà altezza (FWHM) dipende dal solo
dal diametro dell’obiettivo e dalla lunghezza d’onda. Nel caso in cui siano
presenti aberrazioni dinamiche dovute alla turbolenza atmosferica, la PSF
viene deformata sui tempi scala della turbolenza stessa (10−20ms), mostrando peraltro un allargamento della FWHM che determina la diminuzione del
potere risolutivo del sistemo ottico.
Sia dato un oggetto O(x, y) e la PSF del sistema, l’immagine risultante
I(x, y) sarà la convoluzione delle due funzioni:
Z Z
I(x, y) =
O(x, y) · P SF (x, y)dxdy
(1.18)
che nello spazio di Fourier diventa un semplice prodotto convolutivo:
I(u, v) = O(u, v) ∗ P SF (u, v)
(1.19)
tale processo è descritto in figura 1.3.
1.5
Il Seeing diurno
In astronomia, con il termine ”seeing” ci si riferisce all’effetto di aberrazione
delle immagini degli oggetti astronomici dovuto all’atmosfera terrestre Una
9
Figura 1.2: (In alto) PSF teorica del Telescopio Solare THEMIS (90 cm diametro).
(In basso) Simulazione della PSF del Telescopio Solare THEMIS in presenza di alcune
aberrazioni. E’ interessante notare la comparsa di rugosità nella superficie della funzione
causate da ritardi di fase nel fronte d’onda.
Figura 1.3: Processo di convoluzione con la PSF che porta alla formazione dell’immagine
dell’oggetto. La PSF porta ad una degradazione dell’immagine che aumenta all’aumentare
della larghezza dell’FWHM.
10
delle maggiori problematiche dell’ottica adattiva solare risiede nelle turbolenze atmosferiche diurne. L’irraggiamento solare, infatti, è responsabile dell’innesco di forti turbolenze a bassissima quota, con fluttuazioni dell’indice di
rifrazione molto rapide.
Una delle soluzioni a questo problema sta nel collocare l’apertura dello strumento ottico al di sopra dello strato piu’ turbolento, strato che si colloca in
genere tra il suolo e i 10-20 metri di altezza. In figura 1.5 è mostrato il telescopio solare THEMIS (Tenerife, Spagna) con la sua torre che gli permette
oltrepassare lo strato turbolento dovuto al riscaldamento del suolo.
Anche l’intero edificio e la cupola sono stati appositamente studiati per limitare gli effetti delle turbolenze che, a causa dell’irraggiamento, si sviluppano
vicino alle superfici esposte alla luce solare.
Il secondo grande problema delle osservazioni diurne è la turbolenza che
si sviluppa a contatto con le superfici ottiche (“mirror seeing”) e quella che
si sviluppa all’interno dei sistemi ottici stessi. Una delle possibili soluzioni
sta nel realizzare sistemi ottici aperti o a diretto contatto con venti laminari
in grado di portar via le turbolenze superficiali ad alta frequenza. Anche in
questo caso, però, l’uso di alcuni accorgimenti tecnici può solo limitare gli
effetti del seeing e mai cancellarli definitivamente.
La progettazione dei sistemi di compensazione adattiva solare richiede, perciò, particolare cura e interesse verso queste importanti problematiche.
Oltre alle turbolenze a bassa quota e sulle superfici ottiche, esistono le
comuni turbolenze ad alta quota che si manifestano, seppur con minor intensità, anche nelle osservazioni notturne. In Figura 1.4 è riportato un profilo
di velocità del vento in funzione della pressione atmosferica (la quota) per
l’osservatorio del Teide (Tenerife - Isole Canarie, Spagna). Il profilo rappresentato mostra chiaramente un incremento della velocità del vento (fino a
30 m/s) per quote corrispondenti a 200 mbar di pressione. La velocità del
vento in quota rappresenta un elemento fondamentale nella caratterizzazione
delle turbolenze atmosferiche. Profili di questo tipo verranno presi in esame
nel Capitolo 3, quando si parlerà degli algoritmi per la simulazione delle
turbolenze atmosferiche.
1.6
I sensori di fronte d’onda
Il primo passo nella compensazione adattiva è quello dell’analisi e della determinazione del fronte d’onda aberrato.
A tal fine è indispensabile l’impiego di specifici sensori che riescano a deter11
Figura 1.4: Profilo di velocità del vento in funzione della pressione atmosferica per il
sito del Teide (Isole Canarie, Spagna) (Muñoz-Tuñón et al., 2007)
Figura 1.5: Telescopio Solare Themis (Tenerife, Spagna). L’altezza alla quale è posizionato il telescopio permette di ridurre l’incidenza delle turbolenze a bassa quota. Anche la
forma della cupola permette di limitare gli effetti delle turbolenze che si sviluppano nelle
vicinanze delle superfici dell’edificio.
12
Figura 1.6: Schema di funzionamento di un sensore Shack-Hartmann. Un array di lenti,
posto in un piano coniugato alla pupilla, campiona il fronte d’onda in posizioni diverse.
Ogni lente produce un’immagine sul detector. In alto è mostrata la situazione ideale in cui
il fronte d’onda non contiene aberrazioni e le immagini risultanti risultano equispaziate.
In basso è mostrata, invece, la situazione in cui la slope del fronte d’onda provoca lo
spostamento di ogni singola immagine prodotta. Tale spostamento è proporzionale alla
slope locale del fronte d’onda stesso.
minare le fluttuazioni dell’indice di rifrazione che, come visto in precedenza,
risultano associate a variazioni della slope del fronte d’onda dovute a ritardi
di fase.
Uno dei sensori più utilizzati e quindi meglio testati in AO Solare è il sensore
Shack-Hartmann (fig. 1.6). Esso è costituito da un array di lenti poste su
un piano coniugato alla pupilla come mostrato in figura. Ogni lente produce
un’immagine campionando il fronte d’onda in posizioni diverse. In assenza di
aberrazioni (fronte d’onda piatto) le immagini prodotte saranno equispaziate. Al contrario, se il fronte d’onda è corrugato le immagini risulteranno
non equispaziate e ogni spostamento dalla posizione ideale sarà proporzionale
alla slope locale del fronte d’onda.
La presenza di immagini estese e a basso contrasto, quali ad esempio
quelle della granulazione solare, rappresenta una delle maggiori difficoltà nell’impiego di questa tipologia di sensore nell’AO solare.
A differenza del caso notturno, in cui gli spostamenti delle immagini prodotte
risultano dalla misura dello spostamento di una sorgente puntiforme sul piano immagine, nel caso solare la misura degli shift non risulta immediata. Al fine di creare una mappa delle slope del fronte d’onda, è necessario
l’utilizzo di una tecnica di autocorrelazione sulle immagini prodotte dallo
Shack-Hartmann chiamata SAD (Sum of Absolute Differences) (Gelly & the
13
THEMIS team, 2007).
1.7
L’algoritmo SAD
L’algoritmo SAD è un algoritmo ad alta efficienza per la stima del movimento
di un’immagine rispetto ad un’altra immagine di riferimento. Data un’immagine di reference r(i,j) ed un’immagine di test t(i,j) (1 ≤ i ≤ n, 1 ≤ j ≤ n),
l’espressione matematica della SAD è data da:
XX
s(p, q) =
|r(i, j) + t(i + p, j + q)|
(1.20)
La SAD di due immagini identiche è zero nell’origine e negativa altrove.
Quando applicata alle immagini della granulazione solare con spostamenti
random, causati dalla turbolenza atmosferica, la funzione SAD presenterà
un massimo nella posizione dello shift.
La scelta della grandezza della finestra di estrazione della SAD richiede
un compromesso tra massimo campo di ricerca e rapporto segnale-rumore
delle immagini risultanti. La recente indroduzione di istruzioni hardware a
basso livello SSE sui moderni processori, garantisce una elevata facilità di
implementazione della SAD assieme ad elevate prestazioni.
1.8
Gli specchi deformabili
La correzione adattiva richiede l’utilizzo di ottiche deformabili che siano in
grado di compensare la “deformazione” del fronte d’onda aberrato. Facendo
uso delle misure della slope ottenute attraverso lo Shack-Hartmann è possibile, attraverso metodologie di ricostruzione che saranno descritte nel prossimo capitolo, risalire ad una matrice di comandi che applicata allo specchio ne
controlla la sua deformazione. Esistono varie tipologie di specchi deformabili
in base al principio fisico sfruttato per causare la deformazione. I più utilizzati sono: gli specchi bimorfi composti di due lamine di materiali diversi che
si deformano all’applicazione di una tensione ai loro capi, piezoelettrci e a
membrana elettrostatica.
1.9
Schema di controllo e chiusura del loop
di correzione
L’Ottica Adattiva è un sistema in controreazione e come tale deve lavorare in
loop chiuso (Fig. 1.8): il sensore di fronte d’onda misura le deviazioni residue
14
Figura 1.7: Esempio di applicazione della SAD ad una sottopupilla del sensore SH. Il
massimo della superficie corrisponde allo shift relativo dell’immagine dovuto alla variazione
di slope del fronte d’onda in due istanti di tempo successivi. Se il fronte d’onda non fosse
stato aberrato, la SAD avrebbe mostrato un picco al centro della finestra di interesse e
l’intera superficie sarebbe risultata simmetrica.
dalla aplanaticità e conseguentemente invia appropriati comandi allo specchio
deformabile. Il sensore di fronte d’onda viene comunemente posto dopo lo
specchio deformabile. Ciò significa che alla chiusura del loop, il sensore di
fronte d’onda misurerà la fase residua non corretta dallo specchio. Sia data,
infatti, la fase del fronte d’onda aberrato indicata con ϕturbulence (x, y, t). Se
la correzione introdotta dallo specchio deformabile è ϕcorr (x, y, t), la fase
residua si può scrivere come:
ϕres (x, y, t) = ϕturbulence (x, y, t) − ϕcorr (x, y, t)
(1.21)
A loop chiuso, il sistema di AO minimizza la fase residua ϕres (x, y, t), con
risoluzione spaziale e temporale dettata dalla configurazione del WFS. In particolar modo la risoluzione spaziale dipenderà dal numero e dalla disposizione
delle sottopupille, mentre la risoluzione temporale dipenderà fortemente dalla
velocità di acquisizione ed elaborazione delle immagini e di risposta degli attuatori. La fase residua non può essere minimizzata arbitrariamente ma risulta essere soggetta alla qualità della ricostruzione del fronte d’onda. Dato infatti il fronte d’onda aberrato considerato in precedenza, la sua ricostruzione
matematica, operata dal sistema di calcolo (blocco CC in figura), risulterà
15
Figura 1.8: Schema di funzionamento dell’Ottica Adattiva. Il fronte d’onda distorto
incontra lo specchio deformabile e successivamente il sensore che lo analizza. Nelle fasi
iniziali lo specchio deformabile viene mantenuto piatto e il loop è aperto. Dopo il calcolo
della matrice dei comandi, lo specchio viene deformato e il sistema comincia a minimizzare
il segnale di errore entrando in controreazione e chiudendo il loop.
in una stima del fronte d’onda stesso ϕrec (x, y, t). L’errore commesso nella
ricostruzione, comunemente chiamato “fitting-error ” sarà dato da:
ϕf itting−error (x, y, t) = ϕturbulence (x, y, t) − ϕrec (x, y, t)
(1.22)
risulta chiaro allora che:
ϕres (x, y, t) > ϕf itting−error (x, y, t))
(1.23)
dove abbiamo idealmente assunto che ϕres (x, y, t) ≈ ϕcor (x, y, t), che corrisponde al caso in cui lo specchio deformabile riesce a riprodurre perfettamente il fronte d’onda ricostruito.
Anche nel caso in cui il loop di correzione riuscisse a minimizzare perfettamente ϕres (x, y, t), quest’ultimo non potrebbe mai essere nullo ma al
massimo pari al ϕf itting−error (x, y, t) che rappresenta un limite instrinseco del
criterio di ricostruzione adottato. Nel capitolo seguente si analizzeranno in
16
Figura 1.9: Schema del loop di controllo. La fase del fronte d’onda viene campionata
dal sensore di fronte d’onda (WFS) e successivamente ricostruita dal sistema di calcolo
(CC, Control Computer) che elabora la matrice di comandi da inoltrare allo specchio
deformabile (DM)
dettaglio le tecniche di minimizzazione del fitting error basate sull’uso di
metodi PCA (Principal Component Analisys).
Un’ulteriore elemento di disturbo nella correzione AO è rappresentato
dal ritardo temporale. Immaginiamo di ottenere una misurazione del fronte
d’onda al tempo t, la conseguente correzione sarà effettuata solamente al
tempo t + ∆t, dove ∆t rappresenta l’intervallo di tempo di calcolo necessario
per ottenere i coefficienti di deformazione dello specchio. Tale intervallo di
tempo può essere molto lungo se paragonato ai tempi caratteristici della turbolenza atmosferica (10 ms) e ciò rappresenta un forte limite alla qualità
della correzione adattiva
Lo specchio deformabile, quindi, corregge le aberrazioni presenti in atmosfera
misurate ad un istante di tempo precedente a quello di correzione e questo
processo va sotto il nome di “time delay error” (TDE).
Nei capitoli seguenti si analizzeranno alcune tecniche di riduzione di questo
errore basate sulla previsione statistica dei coefficienti di deformazione.
Un altro effetto simile a quello appena descritto è rappresentato dalla variazione del fronte d’onda durante il tempo di integrazione delle immagini del
WFS. La deformazione dello specchio avviene dunque durante l’acquisizione
dell’immagine successiva da parte dello SH, quest’ultima perciò non può essere impiegata per la correzione del fronte d’onda in quanto tiene traccia delle
deformazioni del movimento dello specchio, e questo introduce aberrazioni
che non hanno nulla a che fare con quelle atmosferiche. Tutti questi effetti
concorrono ad incrementare ϕres (x, y, t) rendendolo maggiore del fitting error
17
Figura 1.10: Rappresentazione schematica degli step temporali necessari alla correzione.
La prima immagine viene acquisita e successivamente inviata al sistema di calcolo che
con un certo ritardo elabora i parametri necessari allo specchio per la deformazione. La
deformazione avviene durante l’acquisizione della seconda immagine che quindi non può
essere a sua volta utilizzata per il calcolo dei coefficienti. I tempi mostrati nella sequenza
sono stati stimati per il sistema di correzione previsto per il telescopio solare Themis
(Tenerife - Isole Canarie, Spagna)
in ogni caso. In formule si può scrivere:
ϕres (x, y, t) = ϕf itting−error (x, y, t) + ϕDM error (x, y, t) + ϕT imeDelay
(1.24)
dove ϕDM error (x, y, t) rappresenta l’errore introdotto dallo specchio nel rappresentare la fase ricostruita e ϕT imeDelay (x, y, t) rappresenta l’errore introdotto dal ritardo temporale dovuto al tempo di elaborazione dei segnali provenienti dal WFS.
Il TDE del sistema in retroazione, descritto in figura 1.9, riduce drasticamente le capacità del sistema AO. La riduzione del TDE è un modo
molto efficiente di aumentare l’ampiezza di banda del sistema. Una larghezza di banda troppo piccola o una frequenza di acquisizione del WFS troppo
bassa potrebbero rendere impossibile la minimizzazione dell’errore residuo,
compromettendo la chiusura del loop di retroazione e la stabilità stessa del
sistema.
Nell’ultimo capitolo verrà affrontata in dettaglio la chiusura del loop di correzione e il miglioramento della sua stabilità verificando, attraverso simulazioni, il comportamento delle tecniche di ricostruzione basate su PCA che
18
verranno introdotte nel prossimo capitolo.
Gran parte del lavoro di tesi effettuato si basa sull’analisi e l’implementazione di tecniche e metodologie volte alla minimizzazione dell’errore residuo sia con l’utilizzo di opportune base polinomiali per la ricostruzione, sia
facendo uso di tecniche di previsione statistica delle turbolenze a corto range
temporale.
19
20
Capitolo 2
Caratterizzazione modale del
fronte d’onda
Al fine di ricavare la matrice dei comandi per il controllo della deformazione
dello specchio, è necessario disporre di una ricostruzione del fronte d’onda
aberrato facilmente estraibile dalle misure del WFS. In questo capitolo si
analizzeranno le tecniche di ricostruzione modale del fronte d’onda studiando la qualità della ricostruzione e le tecniche di minimizzazione dell’errore
residuo che, come visto nel capitolo precedente, sono essenziali nel raggiungimento della stabilità del loop di correzione.
Esistono due tecniche di ricostruzione modale. La più usata è basata sull’espansione della fase del fronte d’onda in termini dei polinomi ortonormali di
Zernike. Pur essendo una tecnica di facile implementazione, l’espansione in
polinomi di Zernike richiede un notevole tempo di calcolo che, come abbiamo
già visto, è un elemento essenziale nel corretto funzionamento di un sistema
AO.
La seconda tecnica di ricostruzione è basata sulla PCA (Principal Component
Analisys) che porta a definire un nuovo set di polinomi ortogonali definiti
a loro volta come certe combinazioni degli Zernike. Questi ultimi polinomi sono i polinomi di Karhunen-Loeve che, come dimostreremo in seguito,
garantiscono a parità di calcolo una migliore ricostruzione e un minore errore
residuo.
2.1
Espansione del fronte d’onda in polinomi
di Zernike
Per una apertura circolare senza ostruzione centrale, usando coordinate polari
(r, α), i modi Zernike sono definiti come (Noll, 1976):
21
Figura 2.1: I polinomi di Zernike ordinati secondo lo schema di Noll.
 √

 √2cos(mα) 
n + 1Rnm
2sin(mα)


1(m = 0)
(2.1)
(−1)s (n − s)!
r n−2s
s![(n + m)/2 − s]![(n − m)/2 − s]!
(2.2)
Znm (r, α) =
√
dove:
(n−m)/2
Rnm
=
X
s=0
l’indice n è chiamato grado radiale mentre m è la frequenza azimuthale.
I modi Zernike sono ortonormali su una circonferenza di raggio unitario.
Ogni aberrazione della fase ϕ(r), con r vettore posizione sul dominio della
pupilla, può essere espressa come somma dei modi Zernike:
ϕ(r) =
X
j
22
aj Zj (r)
(2.3)
dove j rappresenta il polinomio j-esimo secondo la usuale nomenclatura
riportata in tabella. I coefficienti aj sono dati da:
Z
aj =
W (r)Zj (r)ϕ(r)dr
(2.4)
dove W (r) è una funzione di peso cosı̀ definita:
W (r) =
1/π(r ≤ 1)
0(r > 1)
(2.5)
Siamo interessati a fluttuazioni casuali del fronte d’onda, i coefficienti dell’espansione saranno perciò anch’essi casuali. E’ possibile dimostrare (Roddier,
2004) che solo i polinomi di Zernike con la stessa frequenza angolare m sono
correlati. Tip-tilt, ad esempio, sono anticorrelati con il termine di coma e il
defocus è anticorrelato all’aberrazione sferica. Le fluttuazioni di fase RMS
del fronte d’onda sono definite come:
2
Z
σ =h
W (r)ϕ2 (r)dri
(2.6)
considerando l’ortogonalità dei polinomi su un cerchio di raggio unitario e
facendo uso della (2.3), la relazione precedente diviene:
σ2 =
X
ha2j i
(2.7)
j
se consideriamo il termine con j = 1 (termine di piston) quest’ultima quantità è infinita. Nel nostro caso siamo interessati solo alle deviazioni del fronte
d’onda da una superficie piatta, possiamo quindi far partire la sommatoria
da j = 2, essendo il termine di piston un semplice offset. I sensori di fronte
d’onda non sono sensibili al termine di piston, questo è il motivo per cui nel
seguito non includeremo questo polinomio nelle analisi e nelle ricostruzioni
del fronte d’onda.
La teoria sopra esposta risulta valida per polinomi definiti su una circonferenza di raggio unitario. L’uso di questi stessi polinomi per descrivere fronti
23
d’onda su pupille anulari è molto rischioso (Hou et al., 2006). All’aumentare
dell’oscuramento centrale, infatti, cresce il cross-talk tra i polinomi e con esso
l’errore nella descrizione delle aberrazioni.
Ricordiamo a tal proposito che i moderni telescopi solari e quelli in fase di
progettazione (EST) sono spesso basati su pupille anulari.
Partendo da un algoritmo per la creazione dei polinomi di Zernike ideato da
Carsten Denker (http://solar.njit.edu/ cdenker/idl/zernike.pro), sono state
effettuate delle modifiche al codice sorgente al fine di includere un parametro
di oscuramento centrale variabile e di includere una procedura di ortogonalizzazione dei polinomi anulari attraveso la tecnica di Gram-Smidth. Attraverso
tale codice, riportato in appendice, è stato infine possibile ottenere un set di
polinomi di Zernike definiti su un dominio anulare adatto alla ricostruzione
modale del fronte d’onda per il telescopio THEMIS.
2.2
Ricostruzione della fase
A partire dal vettore S contenente le M misure degli shift del WFS, è necessario risalire al vettore ϕ contenente N elementi. Il problema può essere
formulato attraverso la relazione lineare:
S = Bϕ
(2.8)
dove B rappresenta la matrice di ricostruzione e M > N, ossia il numero
di incognite (il numero di modi) deve essere inferiore al numero delle misurazioni della slope del fronte d’onda. Più in generale, la ricostruzione della
fase rappresenta quindi un problema di inversione: come stimare la fase a
partire dalle misure del WFS.
Tale processo è, nella pratica, portato a termine mediante la tecnica SVD
(Singular Value Decomposition) la quale fornisce, attraverso la matrice pseudoinversa, l’inversione della relazione (2.8):
ϕ = B −1 S
(2.9)
Nei casi reali un’altra complicazione si aggiunge al quadro appena descritto.
La relazione (2.8) deve in realtà tener conto di un termine aggiuntivo n che
rappresenta un rumore additivo a media zero (Rousset, 1999):
24
S = Bϕ + n
(2.10)
n può essere considerato per semplicità rumore non correlato alle incognite ϕ e con statistica gaussiana.
2.3
PCA: i polinomi di Karhunen-Loeve
I polinomi di Zernike sono un’ottima base per la decomposizione del fronte
d’onda. Molto spesso, però, nei sistemi AO reali si ha a che fare con la minimizzazione dei tempi di calcolo che, come già accennato, sono alla base della
qualità della compensazione. A tal fine è utile definire una nuova base di
polinomi che riesca a garantire tempi di calcolo minimi pur garantendo una
elevata precisione nella ricostruzione. Tale base è rappresentata dai polinomi
di Karhunen-Loeve che hanno anche il grande vantaggio di essere statisticamente indipendenti. Questa rappresentazione modale è definita attraverso la
tecnica PCA.
La PCA è una trasformata su spazio vettoriale capace di ridurre la dimensionalità dei dati in analisi. Tale tecnica si basa sulla diagonalizzazione della
matrice di covarianza e la relativa decomposizione in autovettori fornisce,
attraverso i pesi rappresentati dai rispettivi autovalori, una classificazione
dei modi più rappresentativi. Tali modi, detti di Karhunen-Loeve, possono
essere espressi come combinazione lineare di polinomi di Zernike (Roddier,
2004). Supponiamo di avere un vettore colonna dalle infinite componenti:




A=



a1
a2
a3
.
.
.








(2.11)
La matrice delle covarianze dei coefficienti aj dello sviluppo (2.3) risulta
essere:
25





E(AAT ) = E 



a1
a2
a3
.
.
.








 [a1 a2 a3 ...]






(2.12)
dove E(..) denota una media sull’ensemble. Dal momento che la matrice
delle covarianze è hermitiana, essa può essere diagonalizzata. Ciò significa
che esiste sempre una matrice unitaria U cosı̀ che UE(AAT )UT è una matrice diagonale. I coefficienti bj dell’espansione di Karhunen-Loeve sono le
componenti del vettore colonna B=UA, e la loro matrice di covarianza è:
E(BBT ) = E(UAAT UT ) = UE(AAT )UT
(2.13)
Questa matrice è diagonale e di conseguenza i coefficienti sono statisticamente
indipendenti. E’ utile rappresentare i coefficienti aj dello sviluppo di Zernike
in funzione degli bj attraverso la relazione inversa A = UT B:
aj =
X
ujk bk
(2.14)
k
facendo ancora uso della (2.3) possiamo scrivere la relazione che lega la fase
ai nuovi polinomi di Karhunen-Loeve:
ϕ(r) =
XX
j
k
ujk bk Zj (r) =
X
bk
X
j
k
ujk Zj (r) =
X
k
bk Kk (r)
(2.15)
P
dove i Kk (r) = j ujk Zj (r) sono i modi di Karhunen-Loeve espressi come
combinazione lineare dei polinomi di Zernike.
Dal momento che U è unitaria, il prodotto scalare si conserva e i modi
di Karhunen-Loeve rappresentano una base ortogonale. Ogni funzione di
Karhunen-Loeve può essere riordinata in base al rispettivo autovalore che
“pesa” la capacità del polinomio di contenere l’informazione del fronte d’onda.
In linea di principio è possibile utilizzare un sottoinsieme di polinomi di
Karhunen-Loeve ottenendo un grado di informazione altrimenti ottenibile
solo prendendo in considerazione molti più modi di Zernike. Questa tecnica ci permette quindi di ricostruire il fronte d’onda con un numero molto
limitato di modi, riducendo sensibilmente il tempo di calcolo.
26
2.4
Ricostruzione modale per il telescopio solare THEMIS
Come accennato in precedenza, nella prima fase del lavoro di tesi è stata fatta una analisi comparativa delle due tecniche sopra descritte, facendo uso di
alcuni dati Shack-Hartmann acquisiti presso il telescopio Themis nell’agosto 2006. Questi dati riguardano due modalità di acquisizione dello ShackHartmann, 36 e 22 sottopupille, con una cadenza temporale di circa 14 ms
a loop aperto. Nel seguito chiameremo questi dati “Themis06” per brevità.
Figura 2.2: Setup strumentale del sensore di fronte d’onda del telescopio THEMIS. Le
frecce indicano gli shift relativi dovuti alla variazione di pendenza del fronte d’onda. Come
evidente, l’uso di pupille anulari non permette il campionamento e la ricostruzione nella
zona centrale della fase del fronte d’onda.
L’intero dataset è stato acquisito mediante l’ausilio di uno specchio Tip-Tilt.
Su questi dati è stata applicata la tecnica SAD (Fig. 2.2) vista nel capitolo
precedente e successivamente sono stati stimati i coefficienti di Zernike utilizzando alcuni algoritmi appositamente sviluppati in lnguaggio IDL nel gruppo
di Fisica Solare dell’Università di Roma “Tor Vergata” (DelMoro, 2006). I
coefficienti di Zernike sono stati il punto di partenza per il lavoro di analisi
e di ottimizzazione della ricostruzione basata sulla PCA.
27
Figura 2.3: Fase misurata presso il telescopio THEMIS a loop aperto, ricostruita con 30
polinomi di Zernike
2.4.1
Contributo al Tip-Tilt dei modi di alto ordine
Data una fase di riferimento ricostruita solo attraverso i modi di alto ordine,
si è cercato di stimare il contributo al tip e tilt di tali modi. Tale effetto,
ampiamente descritto in letteratura (Ten Brummelaar, 1995), è stato studiato
attraverso la proiezione della fase citata sui modi di tip e tilt ottenendo una
∆ϕmax ≈ 10−5waves 1 come mostrato in figura.
Questa stima rappresenta anche una stima del cross talk presente tra i
polinomi di alto ordine e il tip-tilt.
2.4.2
Algoritmo IDL per la ricostruzione modale di
Karhunen-Loeve
A partire da un set di polinomi di Zernike definiti e ortonormalizzati su
un dominio anulare, è stato scritto un codice IDL (vedere in appendice
Karhunen−annular.pro) per la costruzione dei relativi polinomi di KarhunenLoeve utilizzando la sequenza Themis06. Tale algoritmo è una rete neurale (Fig. 2.5) nella quale l’apprendimento è caratterizzato dalla media
1
Nel linguaggio comunemente usato nel campo dell’Ottica Adattiva lo sfasamento viene
2π
indicato in “waves“ indicando con esso la quantità
∆δ, dove δ rappresenta il cammino
λ
ottico e λ la lunghezza d’onda.
28
Figura 2.4: Contributo al tip-tilt da parte dei modi di alto ordine.
sull’ensemble di 200 coefficienti di Zernike che rappresentano i coefficienti
mediati temporalmente su una scala di circa 300ms. Nel nostro caso, ogni
polinomio di Karhunen-Loeve è costruito come una combinazione lineare di
200 polinomi di Zernike, ciò significa che utilizzando ad esempio solo il primo modo, la ricostruzione risultante conterrà al suo interno aberrazioni con
fluttuazioni spaziali ad alta frequenza.
L’immagine 2.6 mostra chiaramente come l’utilizzo di soli 10 polinomi di
Karhunen-Loeve riesca a descrivere aberrazioni di ordine più alto rispetto
allo stesso numero di polinomi Zernike. Ciò sottolinea il grande vantaggio
rappresentato dalla ricostruzione di Karhunen-Loeve che, a parità di numero
di modi e quindi di tempo di calcolo, è capace di rappresentare il fronte d’onda molto più accuratamente.
A questo bisogna aggiungere il fatto, già menzionato, che i polinomi di
Karhunen-Loeve, essendo costruiti a partire dalla diagonalizzazione di una
matrice di covarianza, sono statisticamente indipendenti rendendo la proiezione
della fase su ciascun polinomio di gran lunga più efficiente. I polinomi di
Karhunen-Loeve sono stati creati per convenienza a partire dal modo 4 di
Zernike, ciò significa che il tip-tilt non sono stati inclusi nel set. Tale convenienza deriva dal fatto che generalmente questi due ordini vengono corretti
attraverso uno specchio a loro dedicato.
29
Figura 2.5: Schema a blocchi dell’algoritmo per la creazione dei polinomi di Karhunen-
Loeve
Figura 2.6: (Sinistra) Fase di riferimento descritta con 200 polinomi di Zernike. (Centro)
Descrizione della stessa fase attraverso l’utilizzo di 10 polinomi di Zernike. (Destra) Stessa
fase ricostruita utilizzando i primi 10 polinomi di Karhunen-Loeve. E’ utile notare come
la ricostruzione di Karhunen-Loeve sia enormemente più efficiente a parità di numero di
modi utilizzati.
2.5
Fitting error
Nella realtà non è possibile utilizzare un numero illimitato di modi per descrivere un fronte d’onda ed è per questo che qualsiasi ricostruzione rappresenterà
il fronte d’onda reale a meno di un termine, chiamato “fitting error”, che rap30
presenta l’errore commesso nell’utilizzo di un set finito di modi normali.
Il fitting error σf2itting−error è in realtà la somma di tre componenti indipendenti (ming Dai, 1996): il fronte d’onda residuo non compensato ∆ϕN , l’errore
2
dovuto ai cross-talk σr2 e l’errore di misura del WFS σM
:
2
σf2itting−error = ∆ϕN + σr2 + σM
(2.16)
dove N rappresenta il polinomio n-esimo.
Figura 2.7: Fitting error per i polinomi di Zernike (+) e per i polinomi di KarhunenLoeve (*).
Per stimare ∆ϕN si è costruita una fase di riferimento mediante l’utilizzo
di 200 termini di Zernike, successivamente si è studiato l’andamento dell’errore residuo all’umentare del numero di modi utilizzati. L’errore residuo è
stato cosı̀ definito:
Z
∆ϕN =
|ϕN − ϕref erence |2 dr
31
(2.17)
I risultati sono mostrati in figura 2.7. Come è facile notare, i polinomi di
Karhunen-Loeve hanno un fitting error decisamente minore sin dai primi modi utilizzati. Dopo aver impiegato circa 140 modi, il fitting error relativo alla
rappresentazione di Zernike raggiunge quello dei polinomi di Karhunen-Loeve
e non c’è più alcun vantaggio nella scelta tra le due ricostruzioni modali.
2.5.1
L’errore di campionamento
Dalla (2.13) si può vedere che anche l’errore di campionamento e i cross-talk
contribuiscono al fitting error. Sebbene di quest’ultimo non sia possibile dare
una stima accurata, senza testare in dettaglio il comportamento di un eventuale specchio deformabile, l’errore di campionamento è facilmente stimabile
simulando la lettura del fronte d’onda di reference. Dato un certo fronte
d’onda noto, è possibile simulare la misura sovrapponendo su di esso una
maschera numerica di sottopupille e misurando poi gli shift ottenuti dalla
convoluzione della fase con la PSF di ogni lente che compone il sensore.
Le limitazioni maggiori, dettate dal sensore Shack-Hartmann alla ricostruzione,
derivano soprattutto dalla spaziatura delle pupille e dal loro numero. Il sensore Shack-Hartmann, infatti, permette il campionamento del fronte d’onda
con una certa frequenza spaziale dettata dalla spaziatura delle pupille che
compongono l’array (fig 2.2). La procedura che porta alla stima dei coefficienti dei polinomi si basa sull’inversione della matrice delle misurazioni e alla
conseguente soluzione di un sistema di equazioni lineari.
Utilizzare un numero di modi troppo elevato rispetto a quello delle sottopupille (misurazioni) significa avere un sistema di equazioni sottodimensionato. Ciò è proprio quello che è possibile vedere nel grafico 2.8, nel quale il
fitting error non riesce a raggiungere lo zero bensı̀ si attesta ad un valore
di circa 70 waves2 con ben 200 polinomi utilizzati. Tale comportamento è
plausibilmente causato dalla carenza di misure del fronte d’onda dovuta alle
poche sottopupille utilizzate. Come ipotizzanbile, non è possibile ricostruire
il fronte d’onda con un numero arbitrariamente alto di polinomi perché essi,
all’aumentare dell’ordine j, rappresentano fasi con fluttuazioni a frequenze
spaziali sempre più alte, oltre la frequenza di campionamento del WFS.
Se da una parte un comportamento di questo tipo era facilmente immaginabile, ben più interessante è ciò che succede, nelle stesse condizioni, alla
ricostruzione con i polinomi di Karhunen-Loeve. In questo caso, infatti, il
fitting error si attesta rapidamente attorno ad un valore circa costante e l’aggiunta di ulteriori modi nella ricostruzione non porta vistosi miglioramenti.
In definitiva, l’utilizzo di pochi polinomi di Karhuenen-Loeve (i primi 10)
32
Figura 2.8: Fitting error per i polinomi di Zernike (*) e per i polinomi di Karhunen-Loeve
(+) con simulazione di misura dei coefficienti attraverso un sensore Shack-Hartmann con
30 pupille utili. Si noti come, anche in questo caso, i polinomi di Karhunen-Loeve siano
molto più efficienti già a basso ordine.
rappresenta un ottimo compromesso tra minimizzazione dell’errore residuo e
velocità di calcolo.
2.6
Conclusioni
In questo capitolo abbiamo visto come l’uso della tecnica PCA rappresenti un metodo altamente efficace nella descrizione del fronte d’onda aberrato.
Partendo dai polinomi di Zernike, opportunamente definiti su un dominio anulare, si possono costruire i polinomi di Karhunen-Loeve come combinazione
lineare dei primi. Questo ultimo set modale può essere riordinato secondo una
scala gerarchica dettata dal peso di ogni polinomio. Tali pesi sono rappresentati dagli autovalori della matrice di covarianza degli Zernike e rappresentano
la capacità di ogni modo di ricostruire un certo fronte d’onda. Il grande vantaggio dei polinomi di Karhunen-Loeve rispetto agli Zernike è rappresentato
dalla loro indipendenza statistica che è a tutto vantaggio della qualità della
33
descrizione delle aberrazioni. Ogni polinomio di Karhunen-Loeve è definito come combinazione lineare di un certo numero di polinomi di Zernike,
questo numero è scelto opportunamente in base al campionamento spaziale
del WFS.
Questo approccio permette, già con l’utilizzo di pochi polinomi, di rappresentare fronti d’onda altrimenti descrivibili con molte decine di Zernike. Nelle
applicazioni di compensazione adattiva i tempi di correzione rappresentano
un limite molto stringente e da questo punto di vista risulta vantaggioso poter
restringere il set di funzioni utilizzate nella ricostruzione pur mantenendo un
certo livello di accuratezza. La spaziatura delle sottopupille del WFS e la loro
dimensione pongono dei limiti sul campionamento spaziale del fronte d’onda.
A partire dalle misure degli shift delle immagini dello Shack-Hartmann è necessario risalire ai coefficienti di Zernike (o di Karhunen-Loeve) attraverso la
risoluzione di un sistema di equazioni. Nel caso in cui il sistema lineare risulti
sottodimensionato, la ricostruzione darà luogo a degli errori, come discusso
nella sezione 2.6.
Uno dei problemi principali della ricostruzione di un fronte d’onda è la minimizzazione del fitting error. Nella sezione 2.5, si è visto come all’uso dei
polinomi di Karhunen-Loeve sia associato un fitting error sensibilmente minore già dai primi ordini di correzione. Ciò ci suggerisce come i polinomi
di Karhunen-Loeve riescano a minimizzare più efficacemente il fronte d’onda
residuo.
Sebbene i polinomi di Karhunen-Loeve ordinati per covarianza decrescente
costituiscano già un’ottima base per la descrizione del fronte d’onda, un
sensibile miglioramento può essere ottenuto riordinando la base modale per
informazione congiunta.
Usando la teoria dell’informazione, è infatti possibile definire la quantità di
informazione che ogni polinomio può aggiungere alla ricostruzione di una certa fase. Sulla base di questa quantità, chiamata appunto informazione congiunta, può essere fatta una classificazione dei polinomi che più efficacemente
riescono a descrivere la morfologia del fronte d’onda. Tale ottimizzazione porta ad un abbattimento del tempo di calcolo (a parità di errore residuo) dell’
80% rispetto al normale utilizzo dei polinomi di Zernike, garantendo cosı̀ una
resa sorprendentemente elevata anche a bassi ordini di correzione.
La teoria dell’informazione, basata sul concetto di entropia informativa, costituisce quindi un valido strumento per il “denoising” della “Principal Component Analisys”, riuscendo a determinare le componenti modali che, seppur
significative dal punto di vista della potenza, non riescono a massimizzare
l’informazione trasferita dalla sorgente. Molti polinomi di Karhunen-Loeve
hanno infatti un valore dell’informazione congiunta vicino allo zero e possono essere esclusi dalla ricostruzione della fase del fronte d’onda, consen34
tendo un incremento della rapidità di calcolo che rappresenta un requisito
fondamentale nei sistemi AO.
35
36
Capitolo 3
Simulazioni di turbolenza
multi-layer
Le performance di un sistema di compensazione adattiva possono essere studiate mediante l’analisi dell’accuratezza della ricostruzione del fronte d’onda
ottenuta da simulazioni numeriche di turbolenza. Tali simulazioni, sebbene
facciano uso di approssimazioni anche molto forti, sono un buon punto di
partenza per la verifica del comportamento di un sistema AO e di tecniche
matematiche finalizzate all’incremento dell’accuratezza e della stabilità dell’intero processo di correzione.
A tal fine sono stati realizzati due diversi algoritmi IDL (phasescreen.pro e
MTBS.pro in appendice), per simulare il comportamento del sistema di lettura del fronte d’onda alla presenza di turbolenza stratificata, lungo la linea
di vista.
3.1
Turbolenza
La turbolenza è un particolare regime fluidodinamico caratterizzato da processi caotici ad invarianza di scala che determinano rapide fluttuazioni di
pressione e velocità nel fluido stesso.
La turbolenza venne studiata inizialmente da O. Reynolds il quale, alla fine
del XIX secolo, osservò ed investigò scientificamente la transizione tra flussi
laminari e turbolenti. Egli si accorse che un flusso diventava turbolento ogni
qual volta uno specifico parametro, combinazione di velocità U, largezza di
scala L e viscosità cinematica ν, cresceva.
Tale parametro, chiamato appunto “numero di Reynolds”, è definito nel
seguente modo:
37
Re = UL/ν
(3.1)
A bassi numeri di Reynolds (Re < 2300) il flusso è regolare, non appena
Re cresce oltre un limite di circa 4000, il flusso diviene turbolento.
Dalla (3.1) si vede che i flussi laminari possono essere mantenuti tali solo in
presenza di piccole velocità o se essi si sviluppano su piccola scala. Nel 1992
Richardson propose che il processo fosse invariante di scala. Egli immaginò
per la prima volta che “eddies” su grande scala, dovuti a forze esterne agenti
sul flusso, diventassero instabili e si rompessero portando alla formazine di
“eddies” via via più piccoli fino ad arrivare alle scale dominate dalla viscosità
cinematica.
Kolmogorov nel 1941 introdusse il concetto di cascata dissipativa. Egli,
basandosi sull’invarianza di scala, calcolò l’indice della power law associata, supponendo che l’energia fosse iniettata su grande scala, e che fosse poi
trasferita attraverso i processi di frammentazione degli “eddies” su scala sempre più piccola fino ad essere dissipata sulle scale dominate dalla viscosità.
L’ipotesi di Kolmogorov è che la cascata sia essenzialmente unidimensionale,
auto-similare e uniformamente distributa. Il flusso di energia è indipendente
quindi dalla dimensione degli “eddies” alla quale è misurato.
Queste assunzioni sono sufficienti a stabilire, in prima approssimazione, un
indice spettrale e una scala di dissipazione caratteristica (Jimenez, 2004).
Indicando un “eddy” attraverso la sua scala spaziale l e la sua velocità ul ,
l’energia per unità di massa dell’eddy di dimensione l è proporzionale a u2l , e
decade in un tempo dell’ordine di Tl = l/ul , in assenza di dissipazione viscosa.
Il rate con il quale l’energia viene trasferita nel processo di frammentazione
è quindi proporzionale a:
ε∼
u2l
u3
= l
Tl
l
(3.2)
per flussi in equilibrio, questa quantità eguaglia il rate di dissipazione viscosa.
Stiamo assumendo che la statistica del flusso vari molto più lentamente del
rate di energia, in modo tale che esso sia indipendente da l. L’equazione (3.2)
porta direttamente ad una power law per gli incrementi di velocità:
ul ∼ (εl)1/3
38
(3.3)
In questi passaggi abbiamo assunto che il tempo scala di decadimento di un
eddy Tl sia molto minore del tempo scala di dissipazione viscosa l2 /ν. Ciò è
vero fintanto che il numero di Reynolds è alto, ma non appena l’energia fluisce
su scale più piccole questo parametro decresce fino a diventare dell’ordine di
1:
l ∼ η = (ν 3 /ε)1/4
ul l/ν ∼ 1,
(3.4)
questa formula definisce quindi la scala di transizione tra il regime turbolento
e quello dissipativo.
In flussi omogenei, la velocità media tra due punti è zero ma possiamo definire
la sua magnitudine attraverso una funzione di struttura:
S2 (l) = u2l = h[u(x + l) − u(x)]2 i
(3.5)
questa tipologia di relazione appartiene ad una classe ben più ampia di
funzioni, le quali rappresentano i momenti statistici delle funzioni di distribuzione delle differenze di velocità. La funzione di struttura descrive quindi le proprietà statistiche del campo di velocità
La distribuzione di velocità tra “eddies” di diverse dimensioni può essere
convenientemente espressa in termini di spettro di Fourier. Consideriamo un
flusso spazialmente periodico in un box di dimensione L3 , il numero d’onda
associato ad una lunghezza di scala l sarà κ = 2π/l, e le componenti della
velocità, uj , sono espresse in termini di serie o integrali di Fourier del tipo:
u(x) = 2πL−1
X
ũm exp(2πimx/L)
m
+∞
Z
≈
−∞
ũ(κ)exp(iκx)dκ
(3.6)
Per il teorema di Parseval possiamo scrivere:
∗
−1
hu(x)v (x)i = 2πL
39
Z
+∞
−∞
ũ(κ)ṽ ∗ (κ)dκ
(3.7)
Figura 3.1: Rappresentazione schematica della teoria di Kolmogorov attraverso lo spet-
tro energetico spaziale (Falgarone & Passot, 2003). L’energia viene iniettata su grandi scale
spaziali “Drive range” e successivamente trasferita attraverso frammentazione di eddies
verso le scale dove avviene la dissipazione “Dissipation range”. La gamma di scale spaziali
dove avviene il trasferimento di energia verso scale sempre più piccole viene chiamato
“Range inerziale”.
Dal momento che l’energia E ∝ u2l , se u = v è una funzione reale, possiamo interpretare |ũ|2 come lo spettro di energia (Lesieur, 1990). Nel caso
unidimensionale possiamo ad esempio scrivere:
hu21 (x)i
Z
=
+∞
−∞
E11 (κ1 )dκ1
(3.8)
più in generale, nel caso 3-D si ottiene:
Z
hui (x)uj (x)i =
40
+∞
−∞
Eij (κ)dκ
(3.9)
Tutte le funzioni di struttura indrodotte nel capitolo 1 posso quindi essere convenientemente riscritte come decomposizione spettrale nello spazio di
Fourier, in particolar modo è utile scrivere la legge di potenza della cascata
dissipativa:
S2 (l) = u2l = C12 (εl)2/3
(3.10)
in termini di spettro energetico:
−5/3
E11 (κ1 ) = C2 ε2/3 κ1
(3.11)
questa relazione ci suggerisce che la maggior parte dell’energia è a bassi numeri d’onda (grandi scale).
Questa trattazione funziona molto bene su una gamma di scale spaziali intermedie, tra la scala di iniezione e di dissipazione, dette range inerziale.
La Fig.3.1 mostra tre regimi nello spettro energetico che vanno dalle grandi scale spaziali fino alle piccole scale dominate dalla viscosità. Lo spettro
raggiunge un massimo alle scale su cui viene iniettata energia che poi fluisce
attraverso un processo a cascata verso scale via via più piccole. Ad ognuno
di questi tre regimi può essere associato un tempo scala (Lesieur, 2008) che
decresce al diminuire della scala spaziale coinvolta. Grandi eddies, al di sopra della scala di iniziezione, avranno quindi tempi scala lunghi, mentre le
scale diffusive saranno caratterizzate da tempi scala brevi. Il range inerziale,
caratterizzato da processi altamente non lineari, si colloca su tempi scala
intermedi.
3.2
Simulazioni basate su FFT
Esistono principalmente due modi per simulare le turbolenze atmosferiche;
il primo è basato sull’uso dei polinomi di Zernike e del loro spettro a legge
di potenza, il secondo, ben più efficace (Sedmak, 2004), si basa invece sulla
trasformata di Fourier.
La turbolenza atmosferica può essere modellata come un processo random
con scale spaziali e temporali caratteristiche note.
Il grande diametro dei moderni telescopi richiede l’utilizzo di “patch” di turbolenza (detti “phase screen”) di grandi dimensioni. In questo quadro di
riferimento i metodi basati su FFT sono meno soggetti a problematiche quali
ad esempio errori di troncamento, tempo di calcolo e limitate quantità di
41
memoria (Sedmak, 2004).
I metodi di simulazione basati su FFT sono adatti solo a patch quasi monocromatici, a correlazione temporale e single-layer (Welsh, n.d.), ma una situazione più realistica può essere descritta considerando la sovrapposizione di
diversi patch single-layer.
Simulazioni di strati atmosferici di grandi dimensioni (centinaia di metri)
devono essere supportati da matrici di almeno 1024x1024 pixel per ottenere
un campionamento spaziale adeguato. A tale scopo risultano essere molto
validi i metodi introdotti da McGlamery (1976). Il phase screen è ottenuto mediante una trasformata di Fourier inversa del prodotto di una matrice
hermitiana contenente il rumore e la radice quadrata dello spettro della fase.
In applicazioni astronomiche “ground-based ”, lo spettro della fase può essere
approssimato con uno spettro di Von Karman (Welsh, n.d.):
−5/3
Pφ (f ) = 0.00058 r0
(f 2 + 1/L20 )−11/6
(3.12)
dove Pφ è la potenza spettrale, f è la frequenza spaziale, L0 è la “outer scale”
ossia la scala spaziale di iniezione dell’energia e r0 è il parametro di Fried
quasi-monocromatico, definito dalla seguente relazione:
r0 = 0.185 λ
6/5
(cosz)
3/5
Z
∞
0
Cn2 (h)dh
(3.13)
con z angolo zenithale, C la funzione di struttura introdotta nel primo capitolo e h la quota.
Lo spettro di Von Karman degenera in uno spettro di Kolmogorov per L0 che
tende ad infinito. Una volta iniettata dell’energia su una certa scala L0 si assiste ad una cascata dissipativa che coinvolge tutte le scale spaziali inferiori,
fino ad un limite dettato dal dominio delle scale dissipative, certamente inferiori a r0 . Tale limite inferiore di arresto della cascata dissipativa è chiamato
“inner scale” e nel nostro caso può essere trascurato grazie all’elevato indice
della legge di potenza (Sedmak, 2004).
Come si vedrà in seguito, la scala L0 dipende dal layer di turbolenza che
viene preso in considerazione ed in particolar modo dalla quota d’interesse.
I phase screen basati su FFT possono essere simulati attraverso le equazioni
di Johansson & Gavel (1994), modificate opportunamente per correggere gli
effetti di sottocampionamento a bassa ed alta frequenza:
42
Nx /2−1
FT
φFmn
X
=
Ny /2−1
X
Tm′ n′ Hm′ n′ FmF ′FnT′ ·
m′ =−Nx /2 n′ =−Ny /2
· exp[2πi(m′ m/Nx + n′ n/Ny )]
(3.14)
valide sotto le condizioni:
m = −Nx /2, Nx /2 − 1,
n = −Ny /2, Ny /2 − 1
(3.15)
Tm′ n′ = 1,
[m′ /(Nx /2]2 + [n′ /(Ny /2]2 ≤ 1,
(3.16)
Tm′ n′ = 0,
[m′ /(Nx /2]2 + [n′ /(Ny /2]2 > 1,
(3.17)
−5/6
FmF ′FnT′ = 0.1513 (Gx Gy )−1/2 r0
[(m′ /Gx )2 + (n′ /Gy )2 + 1/L0 ]−11/12 (3.18)
FT
dove φFmn
è la fase nello spazio di Fourier, Nx e Ny sono le dimensioni del
supporto in pixel, Tm′ n′ è una maschera binaria antialiasing e Hm′ n′ è una
matrice hermitiana contenente il rumore gaussiano a media 0 e σ = 1.
FmF ′FnT′ è la radice quadrata dello spettro di Von Karman che è nullo per
m′ , n′ = 0, e Gx , Gy è la grandezza del phase screen in metri. In questo caso
si assume che il phase screen sia campionato con passo costante ∆ = ∆x =
Gx /Nx = ∆y = Gy /Ny .
Le simulazioni basate su FFT sono accurate solo per grandi valori del rapporto della dimensione dello screen sulla dimensione del supporto e della
dimensione dello screen su l’outer scale.
In generale, gli effetti del sottocampionamento ad alta frequenza non sono
cosı̀ drammatici a causa del rapido “roll off ” spettrale, mentre gli effetti a
bassa frequenza possono essere ridotti ritagliando una finestra di interesse
nel phase screen.
In figura (3.2) è mostrato un phase screen simulato attraverso il codice IDL
“phasescreen.pro”. Le simulazioni basate su questa tecnica funzionano piuttosto bene su scale spaziali inferiori all’outer scale (Falgarone & Passot,
2003), ossia su quelle scale che denotano il range inerziale della turbolenza
(Fig.3.1). I patch cosı̀ creati devono quindi essere opportunamente apodizzati
a zero su scale maggiori di L0 .
43
Figura 3.2: Phase screen generato attraverso il codice phasescreen.pro. Il phase screen
corrisponde a turbolenza ground layer con r0 = 15 cm, spettro di Von-Karman e outer scale
di 2 metri. Il patch ha dimensioni 4x4 metri con campionamento ∆ = Gx /Nx = 0.0078
3.3
Implementazione multi-layer
I metodi basati su FFT simulano phase screen quasi-monocromatici, single
layer e a correlazione temporale. Nello studio dei sistemi di Ottica Adattiva
si rendono invece necessarie simulazioni con evoluzione temporale.
Inoltre, le simulazioni single layer sono generalmente valide solo per scale
spaziali corrispondenti a turbolenza isotropa, ossia con eddies di dimensioni
inferiori all’outer scale. Di conseguenza, le applicazioni che vedono coinvolti
sistemi ottici di grande apertura devono includere la presenza di anisotropie,
soprattutto se la pupilla è più grande della “ground based outer scale”.
Queste caratteristiche devono essere prese in considerazione attraverso la
combinazione di diversi phase screen single-layer indipendenti (Welsh, n.d.),
ognuno con i propri parametri caratteristici che dipenderanno dalla quota
del layer stesso.
In particolar modo l’evoluzione temporale può essere implementata efficacemente facendo scorrere diversi layer ognuno rispetto all’altro con velocità
intrinseche dettate dalla quota di appartenenza. L’anisotropia stessa può
44
inoltre essere simulata sovrapponendo diversi layer correttamente apodizzati
e spostati l’uno rispetto all’altro.
Figura 3.3: Rappresentazione schematica dell’implementazione multi-layer. Ad ognuno
di essi corrispondono parametri caratterstici differenti scelti in base alla quota di appartenenza. In questo schema x e y rappresentano i vettori che individuano le posizioni sul piano
della pupilla, zk rappresenta la quota dello strato di turbolenza k-esimo l’ungo l’asse ottico
e rs individua la sorgente osservata.
3.3.1
L’algoritmo di simulazione multi-layer
Consideriamo la geometria in figura 3.3. Una serie k di strati di turbolenza
sono disposti parallelamente lungo l’asse z. Il processo random che con45
tribuisce alla fase del fronte d’onda attraverso ciascun layer è Gaussiano,
stazionario e a media zero. Ogni layer, inoltre, contribuisce in modo indipendente alla fase del fronte d’onda risultante.
→
La fase in un punto, individuato dal vettore posizione −
x sul piano (xy)
della pupilla di ingresso del sistema ottico, al tempo t è espressa come
sovrapposizione dei contributi provenienti dai diversi layer:
→
φ(−
x , t) =
K
X
φk
k=1
zk
1−
zs
zk −
−
→
→
x +
r sxy , t
zs
(3.19)
→
dove φk (−
x , t) rappresenta la fase del k-esimo layer, nella posizione trasvera
→
x e al tempo t, zs è l’altezza della sorgente proiettata lungo l’asse z e −
r sxy
−
→
è la proiezione di r s sul piano pupilla.
zk
Nel limite di una sorgente infinitamente distante, il fattore (1 − ) diventa 1
zs
zk −
→
r sxy diventa zk θ, dove θ è l’angolo compreso tra la direzione
e il prodotto
zs
di arrivo dell’onda piana proveniente dalla sorgente e l’asse z.
Invocando l’ipotesi di flusso congelato di Taylor, possiamo esprimere la
variazione temporale del phase screen al tempo t + τ in termini della velocità
trasversa di ogni layer v(zk ), velocità che può essere ricavata dai profili del
tipo riportato in Fig.1.4:
→
φ(−
x , t + τ) =
K
X
k=1
φk
zk
1−
zs
zk −
−
→
→
x +
r sxy − v(zk )τ, t
zs
(3.20)
L’ipotesi di flusso congelato ci permette quindi di interpretare le fluttuazioni spaziali del patch di turbolenza come fluttuazioni temporali.
Il codice IDL da me sviluppato e riportato in appendice (MTBS.pro) implementa questa trattazione simulando la lettura del fronte d’onda derivante da
queste relazioni attraverso una particolare maschera che simula la presenza
di un sensore di fronte d’onda SH. Su ogni apertura della maschera vengono
effettuate delle convoluzioni della PSF terica della sottopupilla con il fronte
d’onda multilayer. Tale procedura, già adottata nell’analisi dati Themis06,
fornisce in output i coefficienti di Zernike della turbolenza simulata.
Il codice MTBS è stato utilizzato per fornire i coefficienti dell’espansioni di
Zernike dell’atmosfera risultante dalla sovrapposizione di due phase screen,
uno posto a dieci metri dal suolo (“Ground Layer G0”) con dimensione
46
10x10m e L0 = 10m e l’altro a 300 metri di quota (layer L1) con dimensioni
300x300m e L0 = 300m. Le velocità trasversali utilizzate sono rispettivamente di 2m/s per G0 e 30m/s per L1, che rispecchiano le tipiche velocità
dei profili misurati.
MTBS implementa, inoltre, la possibilità di scegliere la frequenza di campionamento del sensore SH, tale frequenza è stata scelta in base ai tempi
di coerenza derivanti dai rispettivi parametri di Fried (0.15 cm nel caso in
esame), il tempo tra due campionamenti è quindi pari a 1 ms.
Questa simulazione costituirà l’elemento fondamentale del prossimo capitolo
aprendo la strada allo studio e alla minimizzazione del “Time delay error”
attraverso sequenze dati provenienti da simulazioni.
3.4
Conclusioni
In questo capitolo abbiamo mostrato le procedure di simulazione di patch di
turbolenza atmosferica basate su FFT. Tali simulazioni derivano dalla teoria
di Kolmogorov e seppur accurate nel rappresentare turbolenza al di sotto della scala di iniezione dell’energia non soddisfano i principali prerequisiti utili
allo sviluppo di sistemi di Ottica Adattiva su telescopi di grande diametro.
In particolar modo le simulazioni FFT non possono simulare il comportamento temporale della turbolenza.
Per sopperire a questa mancanza si può far uso dell’approssimazione di Taylor di flusso congelato e considerare la sovrapposizioni di diversi phase screen,
ognuno a diverse quote dal suolo. Ogni phase screen avrà la propria velocità
trasversale e il proprio spettro dettato dalla scala di iniezione di energia e dal
rispettivo parametro di Fried. Attraverso il codice MTBS sviluppato in linguaggio IDL è stata simulata la sovrapposizione di diversi layer di turbolenza
e la relativa lettura da parte di un sensore Shack-Hartmann con un tempo
di campionamento inferiore al tempo di coerenza.
Tale simulazione ha prodotto i coefficienti di Zernike relativi alla turbolenza
in esame con una frequenza di campionamento più elevata della frequenza
di decorrelazione del patch atmosferico. Questi risultati ci consentiranno di
indagare, nel prossimo capitolo, un range temporale inesplorato attraverso
le misure effettuate presso il telescopio THEMIS, per il quale il tempo di
campionamento era maggiore al tempo di coerenza.
47
48
Capitolo 4
Ottimizzazione della sequenza
di correzione attraverso
previsioni di turbolenza basate
su modelli ARMA
Nel primo capitolo abbiamo visto come il TDE 1 sia un elemento fondamentale nella determinazione della frequenza di correzione del sistema AO.
Due elementi concorrono alla determinazione del TDE: in primo luogo i tempi di integrazione che se troppo lunghi (≥ τ0 ) possono includere variazioni
intrinseche della fase del fronte d’onda, fornendo una media temporale del
suo effettivo comportamento, in secondo luogo il fatto che ogni correzione
venga applicata ad un istante di tempo successivo a quello di misura. Oltre
a questo vi è la necessità di escludere dalla sequenza di correzione tutti quei
dati acquisiti durante il movimento dello specchio.
Ci possiamo chiedere a questo punto in che modo sia possibile ridurre questo
effetto.
Per risolvere questo problema negli ultimi anni sono stati proposti ed implementati diversi modelli di previsione a corto termine delle aberrazioni
atmosferiche, basati su algoritmi di previsione lineare o su reti neurali (Jorgenson & Aitken, 1992),(Lloyd-hart & Mcguire, 1995),(Mcguire et al., n.d.),
ottenendo peraltro ottimi risultati in termini di prestazioni (Dessenne et al.,
1999).
Questi modelli hanno come scopo quello di fornire stime delle aberrazioni da
1
In letteratura, in tutti quei casi in cui si fa riferimento ai sistemi AO come sistemi
controreazionati si parla spesso di “Servo Lag error” piuttosto che TDE, intendendo però
lo stesso fenomeno
49
utilizzare al posto delle misure effettuate durante il movimento dello specchio, oppure piccole correzioni che aggiungendosi alle misure reali riescano a
compensare gli effetti di ritardo tra misura e compensazione.
In questo capitolo verrà proposta una tecnica di previsione, ampiamente utilizzata in altri ambiti della ricerca scientifica, che si basa su modelli statistici
ARMA (Auto Regressive Moving Average). A tal fine si renderanno utili le
simulazioni di turbolenza introdotte nel capitolo precedente.
Le misure di fronte d’onda effettuate presso il telescopio Themis, infatti, sono
misure con cadenza temporale maggiore del tempo di coerenza, ed è stato
pertanto impossibile utilizzarle in applicazioni di questo tipo dove è invece
necessario conoscere il fronte d’onda e il suo comportamento su scale temporali molto piccole. Non a caso vedremo che i processi ARMA basano il loro
punto di forza proprio sul comportamento autoregressivo del processo. La
prevedibilità del processo è stata ampiamente verificata sperimentalmente,
mostrando che su scale temporali più piccole del tempo di coerenza vale l’approssimazione di flusso congelato (Jakobsson, 1996), (Roddier et al., 1993).
Il processo statistico associato, inteso come sequenza temporale di coefficienti di Zernike o di Karhunen-Loeve, possiede pertanto caratteristiche autoregressive che introducono un effetto di memoria nella sequenza dati. Tali
proprietà consentono quindi di fare previsioni nel tempo.
4.1
I processi stocastici
Il primo passo nell’analisi delle sequenze temporali consiste nella selezione di
un modello matematico (o una classe di modelli) per i dati in questione. Per
prevedere l’evoluzione temporale dei dati è naturale supporre che ogni osservazione xt sia una realizzazione di una variabile casuale Xt . Ciò ci suggerisce
di modellare i dati come una realizzazione di un processo stocastico, descrivibile attraverso una precisa classe di modelli teorici. Per processo stocastico
si intende una famiglia di variabili casuali {Xt }, definite su un certo spazio
di probabilità Ω.
4.1.1
Stazionarietà
Uno ruolo cruciale nell’analisi di sequenze temporali è ricoperto dai processi
le cui proprietà, o solamente alcune di esse, non variano con il tempo.
Se si desidera fare previsioni, allora dobbiamo necessariamente assumere che
“qualcosa” non varii con il tempo. Nell’estrapolare funzioni deterministiche,
ad esempio, è prassi assumere che la funzione stessa o la sua derivata siano
costanti. L’assumere una derivata costante porta ad una estrapolazione lin50
eare dalla quale segue una previsione.
Nell’analisi delle sequenze temporali il nostro fine sarà quello di prevedere
una serie di risultati che in generale sarà non deterministica, ma conterrà
al suo interno una componente casuale. Se questa componente casuale sarà
stazionaria, si potranno sviluppare potenti tecniche di previsione dei dati.
Più precisamente, una serie temporale {Xt , t = 0, ±1, ...} è stazionaria se alcune delle sue proprietà statistiche coincidono con quelle della serie {Xt+h , t =
0, ±1, ...}, per ogni h intero.
Restringendo la nostra attenzione solo ai momenti di primo e secondo ordine
di {Xt } possiamo definire la stazionarietà in modo più preciso.
Sia {Xt } una serie temporale con E(Xt2 ) < ∞. La media di {Xt } è:
µx (t) = E(Xt )
(4.1)
dove E() rappresenta il valore d’aspettazione. La funzione di covarianza è:
γx (r, s) = Cov(Xr , Xs ) = E[(Xr − µx (r)) − (Xs − µx (s))]
(4.2)
per ogni r e s intero.
Si dice che {Xt } è debolmente stazionaria se (Brockwell & Davis (2002)):
µx (t) indipendente da t
(4.3)
γx (t + h, t) indipendente da t
(4.4)
Possiamo a questo punto definire una funzione di autocovarianza (ACVF
Auto Covariance Function) di {Xt } al tempo ritardato h come:
γx (h) = Cov(Xt , Xt+h )
(4.5)
e una funzione di autocorrelazione (ACF Auto Correlation Function):
ρx (h) ≡
γx (h)
= Cor(Xt , Xt+h )
γx (0)
(4.6)
Nel caso di un processo stazionario si può dimostrare (Brockwell & Davis,
1987) che la ACFV γ(h) gode delle seguenti proprietà:
γ(0) ≥ 0
|γ(h)| ≤ γ(0),
51
(4.7)
∀h ∈ Z
(4.8)
γ(h) = γ(−h),
∀h ∈ Z
(4.9)
la ACVF risulta inoltre definita positiva.
Sebbene le funzioni appena definite siano di grande utilità teorica, in tutti i problemi pratici abbiamo a che fare con sequenze temporali osservate
{x1 , x2 , ..., xn } e non con modelli. Per individuare il modello statistico corrispondente alle nostre osservazioni è molto utile definire la funzione di autocorrelazione campionaria (Sample ACF ).
Siano {x1 , x2 , ..., xn } le osservazioni di una certa serie temporale. La media
campionaria sarà:
1X
x̄ =
xt
n t=1
n
(4.10)
la funzione di autocovarianza campionaria sarà:
n−|h|
1 X
γ̂ =
(xt+|h| − x̄)(xt − x̄),
n t=1
−n < h < n
(4.11)
mentre la funzione di autocorrelazione campionaria sarà:
ρx (h) ≡
γˆx (h)
,
γˆx (0)
−n < h < n
(4.12)
Se consideriamo le osservazioni come una delle possibili realizzazioni della
serie {Xt }, allora la Sample ACF sarà una stima della ACF. Questa stima ci
può suggerire quale sia, dei tanti, il processo stazionario candidato più adatto
a rappresentare la dipendenza funzionale osservata (Box & Jenkins, 1970).
Ogni processo stazionario è infatti caratterizzato da una ACF ben precisa.
Nel seguito ci soffermeremo in dettaglio sui processi ARMA e sulla loro
funzione di autocorrelazione.
4.1.2
Processi ARMA
Una classe di processi stazionari estremamente importante è rappresentata
dai processi ARMA (autoregressive moving average), definiti in termini di
equazioni alle differenze lineari con coefficienti costanti. Essi definiscono una
famiglia parametrica di processi stazionari.
Per ogni funzione di autocovarianza γ tale che limh→∞ γ(h) = 0, e per ogni
intero h > 0, è possibile trovare un processo ARMA con ACVF γx tale
che γx (h) = γ(h). Questo rappresenta uno dei motivi per i quali i processi
ARMA giocano un ruolo determinante nell’analisi delle serie temporali. La
loro struttura lineare, come vedremo, porta ad una teoria della previsione
52
particolarmente semplice ed efficace.
Il processo {Xt } è un processo ARMA(p,q) se è stazionario e se per ogni t,
Xt − φ1 Xt−1 − ... − φp Xt−p = Zt + θ1 Zt−1 + ... + θq Zt−q
(4.13)
dove {Zt } rappresenta un rumore “bianco“ con media 0 e varianza σ 2 , e
φi e θi rappresentano i parametri del modello.
La relazione precedente può essere scritta in forma compatta nel seguente
modo:
φ(B)Xt = θ(B)Zt
t = 0, ±1, ±2, ...,
(4.14)
dove φ e θ sono i polinomi di grado p e q rispettivamente:
φ(z) = 1 − φ1 z − ... − φp z p
(4.15)
θ(z) = 1 + θ1 z + ... + θq z q
(4.16)
e
e B rappresenta l’operatore di shift temporale definito da:
B j Xt = Xt−j
j = 0, ±1, ±2, ...,
(4.17)
I polinomi φ e θ rappresentano rispettivamente la parte autoregressiva e a
media mobile dell’equazione 4.13.
Se φ(z) = 1 allora la 4.13 diventa:
Xt = θ(B)Zt
(4.18)
e il processo è un processo a media mobile di ordine q (MA(q)). E’ abbastanza
evidente che in questo caso l’equazione 4.18 avrà soluzione univoca.
In figura 4.1 è mostrato un esempio di ACF corrispondente ad un modello
MA(1).
Se θ(z) = 1 allora la 4.13 diventa un processo completamente autoregressivo di ordine p (AR(p)), caratterizzato dalla seguente espressione:
53
Figura 4.1: ACF per un modello MA(1) con θ1 = 0.8 (Brockwell & Davis, 2002). Lo
spazio campionario si compone di 200 elementi.
φ(B)Xt = Zt
(4.19)
In questo caso, a differenza del processo puramente a media mobile, l’esistenza e l’unicità della soluzione deve essere appurata mediante un’analisi
accurata.
In figura 4.2 sono mostrate le ACF di tre modelli autoregressivi per tre diversi
valori del parametro θ
Dalla 4.13 emerge chiaramente una prerogativa dei processi AR; i campioni al tempo t hanno ”memoria“ di ciò che è accaduto negli istanti di tempo
precedenti. Questa è una delle proprietà che, come vedremo nel seguito, risulta molto utile nel fornire previsioni basandosi sul campione precedentemente
registrato.
Nella figura 4.3 è mostrato l’andamento della Sample ACF (tratto continuo) stimata attraverso dei dati di esempio (Brockwell & Davis, 2002).
In sovrapposizione è mostrata l’ACF di un modello AR(1) (asterischi) con
θ1 = 0.791, dedotto attraverso una stima parametrica sulla Sample ACF. Si
54
Figura 4.2: ACF per un modello AR(1) per tre diversi valori del parametro θ1 . grafico
(a) θ1 = 0.8, grafico (b) θ1 = 0.3, grafico (c) θ1 = −0.8 (Chatfield (1996)).
noti che il modello di ACF giace al di fuori delle bande di confidenza al 95%
per ritardi che vanno da 2 a 6. Ciò ci suggerisce che un modello AR(1) non è
sufficiente a descrivere il processo caratterizzato dai dati sperimentali, sono
pertanto necessari polinomi autoregressivi di ordine superiore.
Come sottolineato in precedenza, ad ogni processo stocastico corrisponde
una particolare funzione di autocorrelazione. Se ne può concludere, quindi,
che dalla ACF si possano stimare i parametri del modello e da questo ottenere
previsioni dei dati ad istanti di tempo successivi, invertendo la 4.13.
4.2
Previsione di un processo stazionario nel
dominio del tempo
In questo paragrafo analizzeremo il problema della previsione di valori {Xt , t >
n + 1} di un processo stazionario, in termini degli {X1 , ..., Xn }. Tale analisi
è alla base dell’algoritmo di previsione realizzato in linguaggio IDL e finalizzato all’impiego nel campo dell’ottica adattiva solare, riportato più avanti.
L’idea di base è quella di utilizzare le osservazioni effettuate ad un istante di
tempo precedente a n e prevedere, attraverso di esse, il comportamento di
{Xt }.
55
Figura 4.3: ACF sperimentale e modello teorico (asterischi) dedotto attraverso una
stima dei parametri dalla Sample ACF. (Brockwell & Davis, 2002).
Prima di addentrarci nei dettagli di questa tecnica è necessario chiarire il
significato del termine ”previsione“. Supponiamo che {Xt } sia un processo stazionario e che di esso si abbiano n osservazioni Xn . Siamo interessati a trovare una funzione, combinazione lineare degli Xn osservati, che
ci fornisca la migliore previsione m(Xn ) di Xn+h , valore della serie dopo
un intervallo di tempo h. Una definizione matematicamente conveniente
di ”migliore previsione“ è quella di funzione di Xn che minimizzi l’errore
quadratico medio:
E(Xn+h − m(Xn ))2
(4.20)
dove Xn+h è il valore vero della serie dopo un intervallo di tempo h dall’ultima osservazione.
Il nostro scopo è quindi quello di trovare una combinazione lineare degli
Xn , Xn−1, ..., X1 che preveda il valore Xn+h , garantendo la minimizzazione
dell’errore quadratico medio.
56
Nel caso di un processo ARMA(p,q) del tipo 4.13, si può dimostrare
facilmente (Brockwell & Davis, 2002) che la previsione m(Xn+1 ) è esprimibile
nel seguente modo:
m(Xn+1 ) =
n
X
θnj (Xn+1−j − m(Xn+1−j ))
se 1 ≤ n < m
(4.21)
j=1
dove m = max(p, q) e:
m(Xn+1 ) = φ1 Xn +...+φp Xn+1−p +
n
X
θnj (Xn+1−j − m(Xn+1−j )) se 1 ≤ n < m
j=1
(4.22)
Più in generale, il migliore stimatore lineare di una previsione a passo
temporale unitario, m(Xn+1 ), è quindi rappresentato da una relazione lineare
di questo tipo:
m(Xn+1 ) = a0 + a1 Xn + ... + an X1
(4.23)
a questo punto rimane da trovare i coefficienti ai in modo tale da minimizzare
l’errore quadratico medio. Confrontando l’espressione precedente con quella
di un processo 4.13, ci si accorge che i coefficienti ai sono proprio i termini del
polinomio autoregressivo associato. Tali termini possono essere efficacemente
stimati dalla sample ACF.
Per chiarire meglio la procedura di previsione mostriamo l’esempio di un
processo ARMA(1,1) definito dalla seguente relazione:
Xt − φXt−1 = Zt + θZt−1
(4.24)
allora la 4.21 e la 4.22 diventano:
m(Xn+1 ) = φXn + θn1 (Xn − m(Xn ))
57
(4.25)
Un altro caso molto interessante è quello dei processi ARMA(p,1) con
θ1 = 0, applicando ancora una volta le relazioni 4.21 e 4.22 si trova che:
m(Xn+1 ) = φ1 Xn + ... + φp Xn+1−p
(4.26)
L’ultimo esempio è l’essenza della procedura implementata per la previsione
dei coefficienti dei polinomi Zernike riportata nel paragrafo seguente. Tutta
la teoria della previsione si riconduce quindi, come sottolineato, alla stima
dei parametri del processo.
Dal momento che ad ogni processo è associata una specifica ACF, essa risulta
essere un buon punto di partenza per questo scopo.
La stima di questi parametri rappresenta una problematica di notevole complessità che peraltro esula dagli obiettivi di questo lavoro di tesi. In questo
contesto ci limitiamo a dire che la stima dei parametri del modello ARMA è
stata effettuata utilizzando la funzione ”IMSL-ARMA“ delle libreria ”Analisys“ del linguaggio di programmazione IDL (IDL-Reference-Guide (2007)),
la quale effettua la stima utilizzando principalmente due algoritmi: quello
dei momenti per la stima della parte a media mobile e l’algoritmo di YuleWalker per la parte autoregressiva. Per una trattazione dettagliata di entrambi gli algoritmi utilizzati si rimanda alla lettura di trattazioni complete
come: Chatfield (1996), Box & Jenkins (1970) e Brockwell & Davis (2002).
4.3
Previsione dei coefficienti di Zernike
Utilizzando le simulazioni di turbolenza con campionamento temporale di
1 ms mostrate nel capitolo precedente, ed in particolar modo il codice MTBS.pro,
si è cercato di stimare il processo ARMA che caratterizza le serie dei coefficienti di Zernike, per cercare di prevedere il loro andamento nel tempo.
Per far questo si è dapprima trovata la Sample ACF della serie temporale di
ogni serie di coefficienti di Zernike derivanti dalla simulazione e su di esse è
stata fatta una stima dei parametri del processo. In figura 4.4 è mostrato
l’andamento della ACF della serie temporale dei coefficienti del polinomio di
Zernike numero 32 ottenuti dalla simulazione del capitolo precedente.
Il nostro intento è quello di analizzare n osservazioni dei coefficienti di
ogni polinomio di Zernike, calcolare la Sample ACF di ognuna di queste serie, stimare i parametri del processo e prevedere i coefficienti di ogni modo
al tempo n + 1.
Nel caso più semplice in cui si abbia bisogno di una previsione a breve termine (1 o 2 step di campionamento temporale) risulta vantaggioso lavorare
58
Figura 4.4: ACF sperimentale della serie temporale dei coefficienti del polinomio di
Zernike numero 32, calcolata su un set di 100 campioni misurati attraverso il codice di
simulazione MTBS con un passo temporale di 1 ms.
con la parte autoregressiva del modello. Ciò è giustificato dall’andamento
dell’ACF che per piccoli ritardi ha un andamento piuttosto monotono.
Una volta stimati i parametri del modello autoregressivo attraverso la
funzione IDL IMLS-ARMA (figure 4.5 e 4.6) è stata implementata la relazione 4.26 per trovare il valore dei coefficienti all’istante n + 1 ottenendo
i risultati riportati in figura 4.7 (pannello centrale), da confrontarsi con i
valori simulati dei coefficienti (pannello superiore). Nel pannello inferiore è
riportato l’errore residuo.
Questo tipo di previsione si basa su un campione osservativo (100 ms)
che scaturisce da simulazioni di misura a loop aperto, senza la presenza di
uno specchio adattivo. In queste condizioni stiamo osservando e prevedendo
i coefficienti della fase reale del fronte d’onda e non la fase residua non corretta dal sistema AO, come nel caso di loop chiuso.
Alla chiusura del loop di controreazione, come discusso nel capitolo 1, le
aberrazioni di fase misurate si riducono in ampiezza a causa della presenza del nodo di controreazione dello specchio. Ciò che si misura è quindi la
59
Figura 4.5: Modello di ACF teorica ricavato dalla stima parametrica del processo
ARMA(40,1).
Figura 4.6: Scarto quadratico medio tra il modello teorico di ACF ricavato dalla stima
parametrica e la Sample ACF. Il modello teorico approssima molto bene la Sample ACF
mostrando la presenza di errori residui molto piccoli (< 10−7 )
fase residua che lo specchio non riesce a compensare, non più quella effettiva
dovuta alla turbolenza.
E’ ragionevole pensare, di conseguenza, che alla chiusura del loop di con60
Figura 4.7: Coefficienti di Zernike simulati al tempo n + 1 (in alto). (centro) Previsione
dei coefficienti di Zernike al tempo n + 1 effettuata tramite l’utilizzo delle n precedenti
osservazioni e con l’utilizzo di un modello ARMA(40,1). (basso) Scarto assoulto residuo
tra la previsione e il dato osservato.
troreazione anche i coefficienti di Zernike si riducano, riducendo peraltro
l’ampiezza della loro fluttuazione temporale e quindi l’errore della previsione.
Iterando il processo di previsione 4.26 a tempi successivi, si può stimare
l’andamento dell’errore residuo che sarà ovviamente una funzione crescente
del tempo come mostrato in figura 4.8.
questo processo di iterazione si basa esclusivamente sulla conoscenza dell’andamento dei coefficienti in un intervallo temporale di 100ms. La previsione ad istanti di tempo successivi (t > 100ms) risentirà perciò dell’accumularsi dell’errore commesso nella previsione di ogni singolo dato, e il suo
propagarsi nel tempo alle previsioni degli step successivi.
La previsione basata sui modelli ARMA può quindi essere utilizzata per
una stima a breve termine dei coefficienti di Zernike.
61
Figura 4.8: Errore di previsione normalizzato in funzione del tempo per il coefficiente
numero 30. La previsione viene effettuata su un campione di 100 ms di misura. L’errore, di conseguenza, cresce come effetto della perdita di memoria del processo a causa
dell’accumularsi dell’errore su ogni singola stima.
Con riferimento alla sezione 1.9, ricordiamo che nella timeline del processo
di correzione non è possibile considerare i dati acquisiti dal sensore di fronte
d’onda durante il movimento dello specchio. Questo problema di fatto dimezza la frequenza di correzione limitando la qualità della compensazione.
Per sopperire a questo problema si può utilizzare la previsione a breve
termine (10ms) precedentemente descritta al posto delle lacune dovute alla
non applicabilità della correzioni (4.9).
In figura è mostrata la sequenza temporale di correzione, l’immagine numero
2 viene acquisita durante la compensazione da parte dello specchio della fase
analizzata attraverso l’immagine 1. Di conseguenza l’immagine 2 non può
essere utilizzata per effettuare una ulteriore correzione, al suo posto possono
quindi essere utilizzate le previsioni dei coefficienti.
La procedura di correzione si baserebbe a questo punto su una sequenza
62
Istante
[ms]
10
20
30
40
50
60
70
80
90
previsione
t
m(Zi , t) [waves]
Zi (t) [waves]
-0.339234
-0.330861
-0.342804
-0.349827
-0.330762
-0.316258
-0.321010
-0.333579
-0.326996
-0.347135
-0.367973
-0.391349
-0.416221
-0.442045
-0.469443
-0.497641
-0.527122
-0.555256
Tabella 4.1: Esempio previsione iterata su istanti di tempo successivi. Nella colonna di
sinistra sono riportati gli istanti di previsione, in quella centrale il valore dei coefficienti
previsto per ogni istante di tempo e in quella di destra il valore reale. Si noti come la
previsione dei coefficienti si discosti sempre più dal valore reale al crescere del tempo.
Figura 4.9: Sequenza temporale analoga a quella descritta nel primo capitolo con l’ag-
giunta della previsione dei coefficienti a colmare la mancanza di dati privi di bias da inviare
allo specchio deformabile. L’immagine 2 dello Shack-Hartmann viene infatti acquisita
durante il movimento dello specchio e la sua analisi non è quindi utilizzabile.
alternata di dati misurati e previsioni, fornendo correzioni alla stessa frequenza di acquisizione del WFS.
Abbiamo visto come l’errore di previsione cresca col tempo, tale problematica potrebbe essere facilmente risolta mediante l’utilizzo di reti neurali che,
implementate su sistemi di calcolo indipendenti, forniscano una riparametriz63
zazione del processo ARMA con una frequenza temporale dettata dalla soglia
di tolleranza dell’errore di previsione. Il compito di questi algoritmi sarebbe
in sostanza quello di aggiornare la Sample ACF e il suo fit parametrico con
nuove misure, in modo tale da fornire al sistema di compensazione meccanica
una previsione dei coefficienti sempre ottimale.
4.4
Conclusioni
In questo capitolo abbiamo visto come sia possibile utilizzare la teoria dei
processi ARMA per prevedere i coefficienti di Zernike a breve termine. Nella
sequenza di correzione dell’ottica adattiva esistono delle lacune determinate
dalla non applicabilità della correzione causata dal fatto che i parametri
necessari derivano da immagini acquisite durante il movimento dello specchio.
Per colmare queste lacune è stato proposto l’utilizzo di tecniche di previsione,
ampiamente utilizzate in altri ambiti della ricerca, basate su modelli statistici
ARMA. La previsione può essere suddivisa in tre passaggi:
1. Stima della Sample ACF da un set di dati reali
2. Stima dei parametri di un modello ARMA(p,q) a partire da un modello
teorico di ACF che ben approssima quella sperimentale.
3. Utilizzo dei parametri ARMA per la previsione dei coefficienti di Zernike
tramite l’utilizzo dei dati osservati
Questo processo potrebbe essere facilmente implementato su una unità di
calcolo parallelo indipendente, che fornisca al sistema di calcolo principale i
parametri del modello ARMA per la previsione, in modo trasversale e con
frequenza inferiore a quella effettiva di correzione.
Il fondamento dell’intera previsione risiede nella stazionarietà del processo statistico che è alla base del comportamento temporale dei polinomi di
Zernike. Tale stazionarietà è quasi certamente garantita su tempi scala brevi, di qui la necessità di aggiornare frequentemente il modello utilizzato sulla
base di nuove misure di fronte d’onda.
L’intera trattazione esposta in questo capitolo si basa su simulazioni numeriche di turbolenza. Sebbene sia stata provata sperimentalmente la prevedibilità del fronte d’onda e dei relativi coefficienti modali, le tecniche sopra
esposte richiedono, come sviluppo naturale, una verifica su dati sperimentali
sia a loop aperto, sia a loop chiuso, garantendo cosı̀ la possibilità di adattare i
modelli e le procedure alle esigenze dettate dalle problematiche sperimentali.
64
Capitolo 5
Studio della stabilità del
processo di compensazione
attraverso simulazioni di
controreazione
Nel primo capito si è visto che un sistema AO può essere pensato come un
sistema di controreazione multi-variabile. Come i sistemi in controreazione
più comuni, esso è composto da un sensore, il WFS, un sistema di controllo,
il sistema di calcolo CC (Control Computer), e un sistema di compensazione,
lo specchio deformabile.
L’intero sistema è quindi progettato per minimizzare la fase residua non compensata dallo specchio. L’input e l’output di tale sistema di controreazione
sono rispettivamente la fase del fronte d’onda aberrato e il residuo della fase
dopo la correzione del DM.
In questo capitolo useremo la sigla “SLS“ (Servo Loop System) per indicare
il sistema di compensazione adattiva inteso come sistema controreazionato.
La progettazione e l’ottimizzazione di un SLS è un problema molto complesso. Esso coinvolge molti aspetti scientifici e ingegneristici come ad esempio
la turbolenza atmosferica, l’ottica, la meccanica, la teoria dei sistemi servoassistiti e la teoria della descrizione del fronte d’onda.
In questo capitolo ci soffermeremo sull’ultimo aspetto menzionato, quello della ricostruzione del fronte d’onda e del suo contributo all’instabilità del loop
di controreazione, riportando i risultati ottenuti attraverso un codice IDL
appositamente sviluppato per simulare il comportamento del nodo di lettura
e ricostruzione del fronte d’onda del sistema SLS.
65
Figura 5.1: Diagramma rappresentativo di un SLS. La fase del fronte d’onda viene
campionata dal sensore di fronte d’onda (WFS) e successivamente ricostruita dal sistema
di calcolo (CC) che elabora la matrice di comandi da inoltrare allo specchio deformabile
(DM). Ogni blocco dello schema introduce un certo ritardo e un errore che sommandosi
determinano il minimo valore della fase residua consentito.
5.1
Il sistema AO come sistema controreazionato
Nella figura ?? abbiamo introdotto uno schema molto semplificato di SLS
che non tiene conto di alcuni importanti fattori che assumono un ruolo determinante nella stabilità del loop. Ogni blocco del sistema SLS, nella realtà,
introduce infatti un errore e un ritardo i cui effetti devono essere accuratamente verificati.
Lo schema di controllo può quindi essere modificato come mostrato in fig.
5.1. Una volta attivato il sistema di correzione, la fase del fronte d’onda
ϕturbulence (x, y, t) viene sintetizzata e sottratta tramite la deformazione dello
specchio alla fase in ingresso. Di conseguenza, nel loop successivo il WFS si
troverà di fronte la fase residua non corretta dallo specchio ϕres (x, y, t). Il
compito del loop è quello di minimizzare questo residuo in modo tale che la
forma dello specchio compensi al meglio la fase del fronte d’onda incidente.
66
Riprendendo la trattazione del capitolo 1 si può scrivere che:
ϕres (x, y, t) = ϕturbulence (x, y, t) − ϕcorr (x, y, t)
(5.1)
A causa della introduzione di un ritardo e di un errore da parte di ogni blocco
del SLS, la fase riprodotta dallo specchio deformabile non potrà essere identica1 a quella reale appartenente al fronte d’onda nell’istante di correzione, e
di conseguenza non sarà mai possibile annullare la fase residua:
ϕres (x, y, t) > ϕT imeDelay (x, y, t)+
+ϕDM error (x, y, t) + ϕW F Serror (x, y, t) + ϕCCerror (x, y, t)
(5.2)
dove ϕT imeDelay (x, y, t), ϕDM error (x, y, t), ϕW F Serror (x, y, t), ϕCCerror (x, y, t)
rappresentano rispettivamente il contributo alla fase da parte del TDE totale, dell’errore introdotto dallo specchio nel riprodurre il fronte d’onda stimato dal sistema di calcolo, l’errore introdotto dal WFS e quello associato
al sistema di calcolo stesso.
In questo capitolo studieremo in dettaglio il contributo dell’errore e del ritardo temporale del sistema di calcolo e del WFS, trascurando quindi il blocco
DM per il quale si assumerà un comportamento ideale.
Nell’assumere un comportamento ideale per il DM, si suppoone quindi che il
corrispettivo blocco non introduca alcun ritardo temporale, ossia che lo specchio si muova istantaneamente, e che le sue autofunzioni siano esattamente
le stesse utilizzate dal sistema di calcolo per fornire una rappresentazione del
fronte d’onda misurato.
Qualora si conoscessero le autofunzioni intrinseche dello specchio e i suoi
tempi di risposta si potrebbe estendere la trattazione anche al suo contributo.
5.2
Simulazione del loop di controreazione
Per studiare le condizioni di stabilità legate al comportamento del WFS e
del sistema di calcolo è stato realizzato un codice IDL per la simulazione del
loop di correzione. Tale codice si basa sulle misure di fronte d’onda acquisite
sul telescopio THEMIS a loop aperto (senza controreazione)
1
Per fase ”identica“ si intende che i ritardi di fase introdotti dallo specchio riescano a
compensare esattamente quelli del fronte d’onda in ingresso sul sistema AO.
67
Figura 5.2: Diagramma di flusso dell’algoritmo sviluppato per la simulazione del loop di
controreazione. Un fronte d’onda incidente reale viene misurato da un simulatore di WFS
che ricava i relativi coefficienti di Zernike. Tale flusso dati viene processato da un nodo di
calcolo che ricostruisce il fronte d’onda. Il flusso dati viene rinviato al WFS dopo essere
stato sottratto al nuovo fronte d’onda incidente.
68
In fig. 5.2 è mostrato il diagramma di flusso del codice di simulazione.
Il fronte d’onda reale incidente viene misurato attraverso un simulatore di
WFS. Tale WFS ricava i coefficienti di Zernike utilizzati per la ricostruzione
sintetica del fronte d’onda. Questo fronte d’onda viene poi sottratto al nuovo
fronte d’onda incidente appartenente ad un istante di tempo successivo. A tal
proposito ricordiamo che il dataset Themis06 è composto da misure di fronte
d’onda con una cadenza temporale di circa 15ms, la simulazione introduce
quindi lo stesso ritardo temporale tra i vari loop del processo. Questo ritardo
è il TDE introdotto nel primo capitolo. La correzione effettuata in un certo
istante temporale si basa, infatti, su misure di fronte d’onda precedenti che
non descrivono accuratamente il fronte d’onda attuale, introducendo un errore di compensazione. Se il TDE è molto elevato il fronte d’onda riprodotto
dallo specchio è diverso da quello incidente che nel frattempo si è modificato.
Ciò significa che in queste condizioni il sistema AO non riesce a chiudere il
loop di controreazione, non riuscendo a minimizzare la fase residua.
Nella figura 2.8 abbiamo visto come alla ricostruzione sia associato un
errore dipendente dal numero di polinomi utilizzati. Fissato il numero di
sottopupille dello Shack-Hartmann, il problema lineare da risolvere per la
stima dei coefficienti di Zernike è ben posto se il numero di polinomi utilizzati è minore o uguale al numero di sottopupille, in caso contrario il sistema
lineare sarà sottodimensionato e le sue soluzioni non saranno univoche.
5.3
Stabilità del sistema SLS in funzione del
numero di polinomi utilizzati per la ricostruzione
Ciò che ci proponiamo di verificare, in questa sezione, è il grado di instabilità
introdotto dalla ricostruzione del fronte d’onda e dal campionamento del
WFS. Attraverso il codice numerico sono stati simulati i loop di correzione
al variare del numero di polinomi di Zernike utilizzati per la ricostruzione
sintetica del fronte d’onda. Ad ogni loop è stato misurato il grado di piattezza del fronte d’onda residuo che rappresenta l’input del loop successivo,
basandosi sulla seguente definizione di fase residua:
R
ϕloop = R
69
ϕ2res
ϕsupport
(5.3)
dove ϕres è stata introdotta nella 5.1 e ϕsupport rappresenta la funzione
pupilla.
Alla chiusura del loop di controreazione, ϕloop viene minimizzato fino al limite imposto dalla 5.2. In condizioni ideali di compensazione perfetta, ϕloop
tenderebbe rapidamente a 0. Nelle condizioni reali, invece, la convergenza
della fase residua verso il minimo non è affatto garantita, alcune instabilità,
infatti, potrebbero portare il sistema fuori dall’equilibrio rendendolo, non
solo non efficiente, ma addirittura controproducente ai fini della risoluzione
spaziale delle immagini.
Sono state eseguite cinque simulazioni per cinque diversi set di polinomi di
Zernike utilizzati, rispettivamente composti da 20, 30, 40, 50, 100 polinomi.
I risultati sono riportati nei grafici 5.3-5.7.
Figura 5.3: Simulazione di chiusura del loop di controreazione con l’utilizzo di un set
modale composto da 20 polinomi di Zernike.
70
Figura 5.4: Simulazione di chiusura del loop di controreazione con l’utilizzo di un set
modale composto da 30 polinomi di Zernike.
Figura 5.5: Simulazione di chiusura del loop di controreazione con l’utilizzo di un set
modale composto da 40 polinomi di Zernike.
71
Figura 5.6: Simulazione di chiusura del loop di controreazione con l’utilizzo di un set
modale composto da 50 polinomi di Zernike.
Figura 5.7: Simulazione di chiusura del loop di controreazione con l’utilizzo di un set
modale composto da 100 polinomi di Zernike.
72
I grafici mostrano chiaramente come l’utilizzo di un numero di polinomi
maggiore del numero di sottopupille dello Shack-Hartmann (30 nel caso in
esame) porti all’innesco di instabilità tali da impedire la convergenza del loop
di controreazione. Un caso limite è quello riportato in fig. 5.7 in cui l’uso
di 100 polinomi di Zernike porta ad una rapida crescita dell’ampiezza delle
fluttuazioni del fronte d’onda residuo, in totale contrasto con ciò che i sistemi
SLS devono fare.
Il processo di chiusura del loop presenta delle oscillazioni iniziali smorzate
la cui ampiezza cresce col crescere del numero di polinomi utilizzati. Ciò
trova spiegazione nella definizione stessa dei polinomi di Zernike. Essi infatti, come visto nel Capitolo 2, riproducono fluttuazioni di fase di ampiezza
crescente al crescere dell’ordine del polinomio. Nel casi in cui la base di polinomi sia di dimensione inferiore a 30 il loop converge entro 150 passi verso
un valore minimo di circa 0.03 waves, interpretabile attraverso la relazione
5.2 come il contributo del TDE e degli errori di campionamento del WFS e
della ricostruzione modale. Nel caso in cui la base abbia dimensione pari a
40, questo minimo è poco più grande, e nonostante il numero di sottopupille
sia inferiore il sistema SLS rimane stabile.
E’ bene osservare come al crescere del numero di polinomi, le oscillazioni del
sistema non siano più di tipo smorzato ma, al limite, possano determinare
la crescita della fase residua in modo del tutto incontrollato, sbilanciando il
loop di correzione verso l’instabilità.
5.4
Conclusioni
I sistemi AO sono dei sistemi controreazionati il cui compito è quello di minimizzare la fase residua non corretta dallo specchio deformabile.
Partendo dai dati Themis06 si è simulato il comportamento del nodo WFS e
del sistema di calcolo della fase, nella chiusura del loop di controreazione
attraverso un codice IDL. Tali simulazioni, seppur non comprendendo le
problematiche dello specchio deformabile, hanno mostrato risultati molto
interessanti. Al variare del numero di polinomi utilizzati per la ricostruzione,
infatti, i criteri di stabilità sembrano variare, evidenziando la presenza di una
soglia dettata dai parametri del WFS. Utilizzando una base modale con dimensione inferiore o uguale al numero di sottopupille del WFS, il processo di
compensazione adattiva mostra inizialmente delle oscillazioni smorzate per
poi convergere ad un valore costannte della fase residua. Se la base modale
ha dimensioni maggiori del numero di sottopupille, cominciano ad innescarsi
delle oscillazioni che portano il sistema lontano dalla stabilità con risultati
73
addirittura peggiori di quelli che si sarebbero ottenuti senza Ottica Adattiva.
Questi risultati ci spingono a pensare che nella progettazione di un sistema
AO sia necessario analizzare in dettaglio questi fattori, estendendo la trattazione a comprendere gli effetti dello specchio, il cui comportamento a sua
volta contribuisce alle instabilità.
74
Conclusioni generali
Il lavoro di tesi svolto si articola in piú parti autonome, ognuna delle quali
è orientata verso un diverso aspetto della ricostruzione del fronte d’onda per
l’Ottica Adattiva (AO).
Partendo dall’approccio puramente matematico della ricostruzione modale
del fronte d’onda, è stata eseguita una analisi comparativa tra due diverse
tecniche di ricostruzione, rispettivamente l’espansione del fronte d’onda in
polinomi di Zernike e in polinomi di Karhunen-Loeve (K-L).
I polinomi di Karhunen-Loeve sono prodotti attraverso l’uso della tecnica PCA come una combinazione lineare degli Zernike. La base modale di
Karhunen-Loeve, ordinata per covarianza, risulta essere piú efficiente nella
ricostruzione del fronte d’onda di quanto non lo siano i polinomi di Zernike.
L’ordinamento per covarianza non è però l’unico possibile. Nel capitolo 2
abbiamo dimostrato che a partire da una formulazione di entropia congiunta
delle immagini sia possibile attribuire ad ogni funzione K-L un peso relativo
legato al contenuto informativo trasportato da ogni singolo polinomio.
Tale analisi porta ad un riordinamento della base modale K-L basato sull’entropia informativa, che risulta essere di gran lunga piú efficiente perfino degli
stessi polinomi ordinati per covarianza.
L’uso di una efficiente base modale non risolve certamente tutti i problemi
dei sistemi AO.
Molto spesso un sistema AO non riesce a sfruttare tutta la sua larghezza
di banda teorica. Poichè la correzione avviene in un istante di tempo successivo a quello di misura, i sistemi AO presenteranno un errore di ritardo
temporale. Nel capitolo 4, abbiamo mostrato come sia possibile minimizzare
questo errore prevedendo, a breve termine, i coefficienti della ricostruzione
del fronte d’onda. Tale studio è stato basato su simulazioni di turbolenza
multi-layer che ci hanno permesso di ottenere sequenze temporali dei coefficienti di Zernike sui quali applicare modelli statistici ARMA. Attraverso
la stima parametrica di tali modelli è possibile stimare una previsione del
processo.
Nel quinto e ultimo capitolo abbiamo infine studiato la stabilità dei sistemi
75
AO visti come sistemi controreazionati, mostrando come essa possa cambiare
in funzione del numero di polinomi impiegati nella ricostruzione.
Tutto il lavoro di tesi è stato svolto avendo come scopo lo studio della compensazione adattiva su telescopi di grande apertura, caratterizzati da una
pupilla anulare. Questa situazione ben si inserisce nel contesto dello studio
di fattibilità del nuovo telescopio solare EST.
76
Bibliografia
Box, G. E. P., & Jenkins, G. M. 1970. Time Series Analysis, Forecasting,
and Control. Holden-Day, San Francisco.
Brockwell, P. J., & Davis, R. A. 1987. Time Series: Theory and Methods.
Springer-Verlag, New York.
Brockwell, P. J., & Davis, R. A. 2002. Introduction to Time Series and
Forecasting, second edition. Springer-Verlag, New York.
Chatfield, C. 1996. The analysis of time series. An introduction. 5th ed.
Chapman and Hall Texts in Statistical Science Series. London.
DelMoro, D. 2006. Internal report.
Dessenne, C., Madec, P.-Y., & Rousset, G. 1999. Sky implementation of
modal predictive control in adaptive optics. Optics Letters, 24(Mar.),
339–341.
Falgarone, E., & Passot, T. (eds). 2003. Turbulence and Magnetic Fields in
Astrophysics. Lecture Notes in Physics, Berlin Springer Verlag, vol. 614.
FRIED, D. L. 1965. Statistics of a Geometric Representation of Wavefront
Distortion. J. Opt. Soc. Am., 55(11), 1427–1431.
Fried, David L. 1990. Time-delay-induced mean-square error in adaptive
optics. J. Opt. Soc. Am. A, 7(7), 1224.
Gelly, B. F., & the THEMIS team. 2007. THEMIS: the 2006-2012 perspective. and beyond. Memorie della Societa Astronomica Italiana, 78,
13–+.
Hou, Xi, Wu, Fan, Yang, Li, & Chen, Qiang. 2006. Comparison of annular wavefront interpretation with Zernike circle polynomials and annular
polynomials. Appl. Opt., 45(35), 8893–8901.
77
IDL-Reference-Guide. 2007. IDL reference manual. RSI.
Jakobsson, Hans. 1996. Simulations of time series of atmospherically
distorted wave fronts. Appl. Opt., 35(9), 1561–1565.
Jimenez, J. 2004. Lecture notes.
Johansson, E. M., & Gavel, D. T. 1994 (June). Simulation of stellar speckle imaging. Pages 372–383 of: Breckinridge, J. B. (ed), Proc. SPIE
Vol. 2200, p. 372-383, Amplitude and Intensity Spatial Interferometry II,
James B. Breckinridge; Ed. Presented at the Society of Photo-Optical
Instrumentation Engineers (SPIE) Conference, vol. 2200.
Jorgenson, Mark B., & Aitken, George J. M. 1992.
Prediction of
atmospherically induced wave-front degradations. Opt. Lett., 17(7),
466–468.
Lesieur, M. 1990. Turbulence in Fluids. Boston, MA: Kluwer Academic
Publishers. Second Revised Edition.
Lesieur, M. 2008. Turbulence in Fluids. Boston, MA: Springer. Forth Revised
Edition.
Lloyd-hart, Michael, & Mcguire, Patrick. 1995. Spatio-Temporal Prediction for Adaptive Optics Wavefront Reconstructors. Page 95 of: Proc.
European Southern Observatory Conf. on Adaptive Optics.
McGlamery, B. L. 1976. Computer simulation studies of compensation of
turbulence degraded images. Pages 225–233 of: Urbach, J. C. (ed), Proc.
SPIE - Int. Soc. Opt. Eng., Vol. 74, p. 225 - 233. Presented at the Society
of Photo-Optical Instrumentation Engineers (SPIE) Conference, vol. 74.
Mcguire, Patrick C., Rhoadarmer, Troy A., Coy, Hanna, Roger, J., & Angel,
P. Linear Zonal Atmospheric Prediction for Adaptive Optics.
ming Dai, Guang. 1996. Modal wave-front reconstruction with Zernike
polynomials and Karhunen-Loeve functions. J. Opt. Soc. Am. A, 13(6),
1218.
Muñoz-Tuñón, C., Varela, A. M., & Fuensalida, J. J. 2007 (Oct.). Recent
results at the Canarian Observatories. Pages 36–46 of: Revista Mexicana de Astronomia y Astrofisica Conference Series. Revista Mexicana de
Astronomia y Astrofisica, vol. 27, vol. 31.
78
Noll, Robert J. 1976. Zernike polynomials and atmospheric turbulence. J.
Opt. Soc. Am., 66(3), 207–211.
Roddier, F. 2004. Adaptive Optics in Astronomy. Adaptive Optics in Astronomy, Edited by François Roddier, pp. 419. ISBN 0521612144. Cambridge,
UK: Cambridge University Press, November 2004.
Roddier, F., Northcott, M. J., Graves, J. E., McKenna, D. L., & Roddier, D.
1993. One-dimensional spectra of turbulence-induced Zernike aberrations:
time-delay and isoplanicity error in partial adaptive compensation. J. Opt.
Soc. Am. A, 10(5), 957.
Rousset, G. 1999. Wave-front sensors. Adaptive Optics in Astronomy. Pages
91–+.
Sedmak, Giorgio. 2004. Implementation of Fast-Fourier-Transform-Based
Simulations of Extra-Large Atmospheric Phase and Scintillation Screens.
Appl. Opt., 43(23), 4527–4538.
Ten Brummelaar, T. A. 1995. The contribution of high order Zernike modes
to wavefront tilt. Optics Communications, 115(Apr.), 417–424.
Welsh, Byron M. A Fourier Series Based Atmospheric Phase Screen
Generator for Simulating Nonisoplanatic Geometries and Temporal
Evolution.
79
80
Appendice sorgenti IDL
pro i n f o r m a t i o n
DEVICE, DECOMPOSED = 0
loadct , 0
; r e s t o r e , ’ / media / d i s k / z e r n i k e 2 0 0 G . sav ’ , / v e r b o s e
; h i s t=h i s t o g r a m ( z e r ( ∗ , ∗ , 8 ) )
restore ,
restore ,
’ / media / d i s k / p h a s e r e f . s a v ’
’ / media / d i s k / K L 2 0 0 f l o a t . s a v ’
; f o r i =0 ,1043 do b e g i n
;
f o r j =0 ,1043 do b e g i n
;
i f p h a s e r e f ( i , j ) l t −10 t he n p h a s e r e f ( i , j )=0 e l s e p h a s e r e f ( i , j )=1
;
endfor
; endfor
p h a s e r e f=p h a s e r e f −mean ( p h a s e r e f )
r e c= f l t a r r ( 2 0 0 )
en tr o p y 2= f l t a r r ( 2 0 0 )
j o i n t E= f l t a r r ( 2 0 0 )
window , 0 , x s i z e =600 , y s i z e =600
tvscl , congrid ( phase ref ,600 ,600)
h i s t=h i s to g r a m ( p h a s e r e f )
h i s t=h i s t / t o t a l ( h i s t )
window , 1 , t i t l e = ’ h i s to g r a m ’
plot , h i s t
en tr o p y =−( t o t a l ( h i s t ∗ a l o g ( h i s t ) ) )
print , en tr o p y
window , 2 , t i t l e = ’ en tr o p y ’
plot , −h i s t ∗ ( a l o g ( h i s t ) )
f o r r =2 ,199 do begin
r e c ( r )= t o t a l ( ( ( karhunen ( ∗ , ∗ , r ) ) ∗ p h a s e r e f )/ t o t a l ( karhunen ( ∗ , ∗ , 0 ) ) )
h i=h i s to g r a m ( r e c ( r ) ∗ karhunen ( ∗ , ∗ , r ) )
h i=h i / t o t a l ( h i )
en tr o p y 2 ( r )=−( t o t a l ( h i ∗ a l o g ( h i ) ) )
j o i n t P=abs ( hi −h i s t )/ h i s t
81
j o i n t E ( r)=− t o t a l ( h i ∗ h i s t ∗ ( a l o g ( j o i n t P ) ) )
j o i n t P=j o i n t P ∗ 0 .
h i=h i ∗ 0 .
endfor
m u t u a l i n f o r m a t i o n=en tr o p y 2−j o i n t E ∗ 0 .
m u t u a l i n f o r m a t i o n ( 0 : 3 ) = en tr o p y 2−j o i n t E ( 0 : 3 )
; m u t u a l i n f o r m a t i o n (4:199)= r e v e r s e ( e n t r opy 2−j o i n t E ( 4 : 1 9 9 ) )
window , 4 , x s i z e =900
plot , m u t u a l i n f o r m a t i o n
loadct ,13
o p l o t , m u t u a l i n f o r m a t i o n , psym=5
stop
o r d e r =( s o r t ( m u t u a l i n f o r m a t i o n ( 4 : 1 9 9 ) ) )
k l=karhunen ∗ 0 .
mi=m u t u a l i n f o r m a t i o n ∗ 0 .
f o r j =0 ,195 do begin
k l ( ∗ , ∗ , j +4)=karhunen ( ∗ , ∗ , o r d e r ( j )+4)
mi ( j )= m u t u a l i n f o r m a t i o n ( o r d e r ( j )+4)
endfor
k l ( ∗ , ∗ , 0 : 3 ) = karhunen ( ∗ , ∗ , 0 : 3 )
stop
m u t u a l i n f o r m a t i o n ( 4 : 1 9 9 ) = mi ( 0 : 1 9 5 )
loadct , 0
window , 4 , x s i z e =900
plot , m u t u a l i n f o r m a t i o n
loadct ,13
o p l o t , m u t u a l i n f o r m a t i o n , psym=5
; save , f i l e n a m e =’/media / d i s k / karhunen MI . sav ’ , k l
stop
end
pro c l o s e d l o o p
r e s t o r e , ’ / media/ d i s k −1/ z e r n i k e 2 0 0 G . s a v ’ , / v e r b o s e
; r e s t o r e , ’ / media / d i s k / t o t z v e c t t s 2 . sav ’ , / v e r b o s e
; t o t z v e c t t s (ZE , image )
82
r e s t o r e , ’ / media / d i s k −1/ t o t z v e c t 2 0 0 . s a v ’ , / v e r b o s e
t o t z v e c t t s=t o t z v e c t
m a x z e r n i k e =100
d e l a y =3
p h a s e r e a l=z e r ( ∗ , ∗ , 0 ) ∗ 0 .
p h a s e0=p h a s e r e a l ∗ 0 .
e r r o r= f l t a r r ( 5 0 0 , 8 )
max=10
k imm=1
n imm=200
Z d i m en s i o n =200
Z v e c t= f l t a r r ( Z d i m en s i o n )
; number o f subap i n t h e SH f i r s t an n u l u s
N 1=12
; number o f subap i n t h e SH s e c on d an n u l u s
N 2=18
; t o t a l number o f s u b a p e r t u r e i n t h e SH ar r ay
N Subap=N 1+N 2
Xm=60
; number o f columns i n each s u b a p e r t u r e
Ym=60
; number o f rows i n each s u b a p e r t u r e
Dim X=712
; number o f rows o f t h e SH ar r ay
Dim Y=632
Centr POS=i n t a r r ( 2 , N Subap ) ; s u b a p e r t u r e c e n t r e p o s i t i o n s i n t h e SH ar r ay
Centr PUP=i n t a r r ( 2 , N Subap ) ; s u b a p e r t u r e c e n t r e p o s i t i o n s i n t h e p u p i l ar r ay
Xb=64
; number o f columns i n s u b a p e r t u r e u s e d f o r c o r r e l a t i o n
Yb=64
; number o f rows i n s u b a p e r t u r e u s e d f o r c o r r e l a t i o n
x c o r r =16
; number o f columns i n c o r r e l a t i o n ar r ay
y c o r r =16
; number o f rows i n c o r r e l a t i o n ar r ay
MAX R=4
;####################################################
; ∗∗∗∗∗∗∗∗∗∗ ∗ ∗ ∗∗ ∗ ∗∗ ∗ ∗ ∗∗ ∗ ∗ ∗∗ ∗ ∗∗ ∗ ∗ ∗∗ ∗ ∗ ∗∗ ∗ ∗∗ ∗ ∗ +
telescope D = 0.90
; 0. 475
; T e l e s c o p e a p e r t u r e w i d t h (m)
;0.0
; T e l e s c o p e c e n t r a l obscurament w i d t h (m)
t e l e s c o p e O = 0 . 4 0 ; 0. 396
lambda
= 5 0 0 0 . e −10
; 4305. e−10
; 4 5 0 7 . e−10
; Working w a v e l e n g t h (m)
pix
=
0.215
; P i x e l −s i z e en a r c s e c ( 0 . 0 8 2 9 3 7 a t 450 nm) ( 0 . 0 8 2 8 5 2 a t 430 nm)
s i z =500 ; 475
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
TOT Z VECT= f l t a r r ( Z d i m en s i o n , max )
SQR Z VECT= f l t a r r ( Z d i m en s i o n , max )
T o t v e c t= f l t a r r (2 ∗ N subap , max)
print , s y s t i m e ( 0 )
r e s t o r e , ’ / media / d i s k −1/G200 . s a v ’ , / v e r b o s e
G0=G( 0 : Z d i m en s i o n −1 ,∗)
G=G0
svdc ,G,W, U,V
83
y1=7+ f i x (ym / 2 . ) ; Cen x −2∗Xm−2∗ d i s t /2
y2=99+ f i x (ym / 2 . ) ; Cen x −1∗Xm−1∗ d i s t /2
y3=192+ f i x (ym / 2 . ) ; Cen x
y4=284+ f i x (ym / 2 . ) ; Cen x+1∗Xm+1∗ d i s t /2
y5=376+ f i x (ym / 2 . ) ; Cen x+2∗Xm+2∗ d i s t /2
y6=468+ f i x (ym / 2 . ) ; Cen x+2∗Xm+2∗ d i s t /2
y7=559+ f i x (ym / 2 . ) ; Cen x+2∗Xm+2∗ d i s t /2
x1=10+ f i x (xm / 2 . ) ; Cen y −2∗Ym−2∗ d i s t
x2=62+ f i x (xm / 2 . ) ; Cen y −1.5∗Ym−1∗ d i s t
x3=115+ f i x (xm / 2 . ) ; Cen y −1∗Ym−1∗ d i s t
x4=168+ f i x (xm / 2 . ) ; Cen y −0.5∗Ym−0.5∗ d i s t
x5=222+ f i x (xm / 2 . ) ; Cen y
x6=274+ f i x (xm / 2 . ) ; Cen y +0.5∗Ym+0.5∗ d i s t
x7=327+ f i x (xm / 2 . ) ; Cen y+1∗Ym+1∗ d i s t
x8=380+ f i x (xm / 2 . ) ; Cen y +1.5∗Ym+1∗ d i s t
x9=434+ f i x (xm / 2 . ) ; Cen y+2∗Ym+2∗ d i s t
x10=486+ f i x (xm / 2 . ) ; Cen y+2∗Ym+2∗ d i s t
x11=539+ f i x (xm / 2 . ) ; Cen y+2∗Ym+2∗ d i s t
x12=592+ f i x (xm / 2 . ) ; Cen y+2∗Ym+2∗ d i s t
x13=646+ f i x (xm / 2 . ) ; Cen y+2∗Ym+2∗ d i s t
; ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 FIRST ANNULUS
Centr POS ( ∗ , 1 8 ) = [ x7 , y6 ]
Centr POS ( ∗ , 1 9 ) = [ x9 , y6 ]
Centr POS ( ∗ , 2 0 ) = [ x10 , y5 ]
Centr POS ( ∗ , 2 1 ) = [ x11 , y4 ]
Centr POS ( ∗ , 2 2 ) = [ x10 , y3 ]
Centr POS ( ∗ , 2 3 ) = [ x9 , y2 ]
Centr POS ( ∗ , 2 4 ) = [ x7 , y2 ]
Centr POS ( ∗ , 2 5 ) = [ x5 , y2 ]
Centr POS ( ∗ , 2 6 ) = [ x4 , y3 ]
Centr POS ( ∗ , 2 7 ) = [ x3 , y4 ]
Centr POS ( ∗ , 2 8 ) = [ x4 , y5 ]
Centr POS ( ∗ , 2 9 ) = [ x5 , y6 ]
; ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 SECOND ANNULUS
Centr POS ( ∗ , 0 ) = [ x8 , y7 ]
Centr POS ( ∗ , 1 ) = [ x10 , y7 ]
Centr POS ( ∗ , 2 ) = [ x11 , y6 ]
Centr POS ( ∗ , 3 ) = [ x12 , y5 ]
Centr POS ( ∗ , 4 ) = [ x13 , y4 ]
Centr POS ( ∗ , 5 ) = [ x12 , y3 ]
Centr POS ( ∗ , 6 ) = [ x11 , y2 ]
Centr POS ( ∗ , 7 ) = [ x10 , y1 ]
Centr POS ( ∗ , 8 ) = [ x8 , y1 ]
Centr POS ( ∗ , 9 ) = [ x6 , y1 ]
Centr POS ( ∗ , 1 0 ) = [ x4 , y1 ]
Centr POS ( ∗ , 1 1 ) = [ x3 , y2 ]
Centr POS ( ∗ , 1 2 ) = [ x2 , y3 ]
Centr POS ( ∗ , 1 3 ) = [ x1 , y4 ]
Centr POS ( ∗ , 1 4 ) = [ x2 , y5 ]
Centr POS ( ∗ , 1 5 ) = [ x3 , y6 ]
Centr POS ( ∗ , 1 6 ) = [ x4 , y7 ]
Centr POS ( ∗ , 1 7 ) = [ x6 , y7 ]
TOT Z VECT=TOT Z VECT ∗ 0 .
p h a s e r e a d=p h a s e r e a l ∗ 0 .
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
; ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ SIMULAZIONE LOOP ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
84
f o r image =2 ,20 , d e l a y do begin
p h a s e r e a l=p h a s e r e a l ∗ 0 .
f o r i =0 , m a x z e r n i k e do begin
p h a s e r e a l (∗ , ∗ )= p h a s e r e a l (∗ , ∗ )+ Zer ( ∗ , ∗ , i )∗ t o t z v e c t t s ( i , image )
p h a s e r e a d (∗ , ∗ )= p h a s e r e a d (∗ , ∗ )+ Zer ( ∗ , ∗ , i )∗ t o t z v e c t ( i , 0 )
endfor
shade surf , phase real
mm=0
; f o r mm=0,1 do b e g i n
p h a s e r e f=p h a s e r e a l −p h a s e r e a d
p h a s e r e a d=p h a s e r e a d ∗ 0 .
s u p p o r t=z e r ( ∗ , ∗ , 0 )
REFERENCE= f l t a r r ( Dim X , Dim Y )
imm=p h a s e r e f
For i s u b a p =0 , N subap −1 do begin
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
xpos=Centr POS ( 0 , i s u b a p )− f i x (xm / 2 . )
ypos=Centr POS ( 1 , i s u b a p )− f i x (ym / 2 . )
sub imm=imm( xpos : xpos+xm−1 , ypos : ypos+ym−1)
REFERENCE( xpos : xpos+xm−1 , ypos : ypos+ym−1) =sub imm
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
endfor
vx=i n t a r r (xm)
vy=i n t a r r (ym)
For x=0 ,xm−1 do vx ( x)=x−xm/2
For y=0 ,ym−1 do vy ( y)=y−ym/2
vx= s h i f t ( vx , xm/ 2 )
vy= s h i f t ( vy , ym/ 2 )
Mx=vx#r e p l i c a t e ( 1 ,ym)
My=r e p l i c a t e ( 1 ,xm)#vy
For x=0 ,xm−1 do vx ( x)=x
For y=0 ,ym−1 do vy ( y)=y
Px=vx#r e p l i c a t e ( 1 ,ym)
Py=r e p l i c a t e ( 1 ,xm)#vy
v s h i f t= f l t a r r ( 2 , N Subap )
; ˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆ
n u c u t o f f=t e l e s c o p e D / lambda
; Fr e c u e n c i a de c o r t e en radˆ−1
pixrad = pix ∗! pi /180./3600.
; P i x e l −s i z e en r a d i a n .
d e l t a n u =1./ s i z / p i x r a d
; I n t e r v a l o de muestreo en radˆ−1
; Radio de l a p u p i l a en p i x e l s
r p u p i l=n u c u t o f f / d e l t a n u / 2 .
85
p i x c u t =2∗ r p u p i l
s i z =2.ˆ f i x ( a l o g 1 0 ( p i x c u t )/ a l o g 1 0 ( 2 . ) +1) +20 ; l a prima m a t r i c e con dim 2ˆn che c o n t i e n e l a
s h i f t x= f l t a r r ( N Subap )
s h i f t y= f l t a r r ( N Subap )
ARRAY= f l t a r r ( Dim x , Dim y )
PSF ARRAY= f l t a r r ( Dim x , Dim y )
PSF ref ARRAY= f l t a r r ( Dim x , Dim y )
D v ect= f l t a r r (2 ∗ N subap )
; rd=radius aper ANN ( s i z /2 , p i x c u t , t e l e s c o p e O / t e l e s c o p e D , t s u p p o r t )
; ˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆ;ˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆ
;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
Cen=f i x ( s i z / 2 . )
e x a g s i d e= 84 ; 7 9 ; 1 0 9 ; 2 0 0
e x a g a p t=e x a g s i d e ∗ c o s ( ! PI / 6 . )
; ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12
y1=Cen −4.5∗ e x a g s i d e
y2=Cen−3∗ e x a g s i d e
y3=Cen −1.5∗ e x a g s i d e
y4=Cen
y5=Cen +1.5∗ e x a g s i d e
y6=Cen+3∗ e x a g s i d e
y7=Cen +4.5∗ e x a g s i d e
; ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12
x1=Cen−6∗ e x a g a p t
x2=Cen−5∗ e x a g a p t
x3=Cen−4∗ e x a g a p t
x4=Cen−3∗ e x a g a p t
x5=Cen−2∗ e x a g a p t
x6=Cen−1∗ e x a g a p t
x7=Cen
x8=Cen+1∗ e x a g a p t
x9=Cen+2∗ e x a g a p t
x10=Cen+3∗ e x a g a p t
x11=Cen+4∗ e x a g a p t
x12=Cen+5∗ e x a g a p t
x13=Cen+6∗ e x a g a p t
; ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12
;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
;++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 FIRST ANNULUS
Centr PUP ( ∗ , 1 8 ) = [ x7 , y6 ]
Centr PUP ( ∗ , 1 9 ) = [ x9 , y6 ]
Centr PUP ( ∗ , 2 0 ) = [ x10 , y5 ]
Centr PUP ( ∗ , 2 1 ) = [ x11 , y4 ]
Centr PUP ( ∗ , 2 2 ) = [ x10 , y3 ]
Centr PUP ( ∗ , 2 3 ) = [ x9 , y2 ]
Centr PUP ( ∗ , 2 4 ) = [ x7 , y2 ]
Centr PUP ( ∗ , 2 5 ) = [ x5 , y2 ]
Centr PUP ( ∗ , 2 6 ) = [ x4 , y3 ]
Centr PUP ( ∗ , 2 7 ) = [ x3 , y4 ]
Centr PUP ( ∗ , 2 8 ) = [ x4 , y5 ]
Centr PUP ( ∗ , 2 9 ) = [ x5 , y6 ]
; ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 SECOND ANNULUS
Centr PUP ( ∗ , 0 ) = [ x8 , y7 ]
Centr PUP ( ∗ , 1 ) = [ x10 , y7 ]
Centr PUP ( ∗ , 2 ) = [ x11 , y6 ]
Centr PUP ( ∗ , 3 ) = [ x12 , y5 ]
Centr PUP ( ∗ , 4 ) = [ x13 , y4 ]
Centr PUP ( ∗ , 5 ) = [ x12 , y3 ]
86
Centr PUP ( ∗ , 6 ) = [ x11 , y2 ]
Centr PUP ( ∗ , 7 ) = [ x10 , y1 ]
Centr PUP ( ∗ , 8 ) = [ x8 , y1 ]
Centr PUP ( ∗ , 9 ) = [ x6 , y1 ]
Centr PUP ( ∗ , 1 0 ) = [ x4 , y1 ]
Centr PUP ( ∗ , 1 1 ) = [ x3 , y2 ]
Centr PUP ( ∗ , 1 2 ) = [ x2 , y3 ]
Centr PUP ( ∗ , 1 3 ) = [ x1 , y4 ]
Centr PUP ( ∗ , 1 4 ) = [ x2 , y5 ]
Centr PUP ( ∗ , 1 5 ) = [ x3 , y6 ]
Centr PUP ( ∗ , 1 6 ) = [ x4 , y7 ]
Centr PUP ( ∗ , 1 7 ) = [ x6 , y7 ]
; ı̈¿ 21 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21
;++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Exag mask=i n t a r r ( e x a g s i d e ∗2+1 , e x a g s i d e ∗2+1)
window , 3 , x s i z e=e x a g s i d e ∗2+1 , y s i z e=e x a g s i d e ∗2+1 , t i t l e = ’ Exag mask ’
t v s c l , Exag mask
xpos=e x a g s i d e
ypos=e x a g s i d e
y1=Ypos−e x a g s i d e
y2=Ypos−e x a g s i d e /2
y3=Ypos+e x a g s i d e /2
y4=Ypos+e x a g s i d e
x1=Xpos−e x a g a p t
x2=Xpos
x3=Xpos+e x a g a p t
p l o t s , [ X2 , X3 ] , [ Y1 , Y2 ] , / d e v i c e , c o l o r = [ 2 5 5 , 2 5 5 , 2 5 5 ] , Thick=2
p l o t s , [ X3 , X3 ] , [ Y2 , Y3 ] , / d e v i c e , c o l o r = [ 2 5 5 , 2 5 5 , 2 5 5 ] , Thick=2
p l o t s , [ X3 , X2 ] , [ Y3 , Y4 ] , / d e v i c e , c o l o r = [ 2 5 5 , 2 5 5 , 2 5 5 ] , Thick=2
p l o t s , [ X2 , X1 ] , [ Y4 , Y3 ] , / d e v i c e , c o l o r = [ 2 5 5 , 2 5 5 , 2 5 5 ] , Thick=2
p l o t s , [ X1 , X1 ] , [ Y3 , Y2 ] , / d e v i c e , c o l o r = [ 2 5 5 , 2 5 5 , 2 5 5 ] , Thick=2
p l o t s , [ X1 , X2 ] , [ Y2 , Y1 ] , / d e v i c e , c o l o r = [ 2 5 5 , 2 5 5 , 2 5 5 ] , Thick=2
b=LABEL REGION(255 − t v r d ( ) , /EIGHT)
Exag mask=b eq 3
w d el ete , 3
; ˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆ;ˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆ
;FOR Z J =2, Z dimension −1 DO BEGIN
Z v e c t=Z v e c t ∗0
REFERENCE= f l t a r r ( Dim X , Dim Y )
imm=p h a s e r e f
For i s u b a p =0 , N subap −1 do begin
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
xpos=Centr POS ( 0 , i s u b a p )− f i x (xm / 2 . )
ypos=Centr POS ( 1 , i s u b a p )− f i x (ym / 2 . )
sub imm=imm( xpos : xpos+xm−1 , ypos : ypos+ym−1)
REFERENCE( xpos : xpos+xm−1 , ypos : ypos+ym−1) =sub imm
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
endfor
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
d e v i c e , decomposed=0
l o a d c t , 3 , / SILENT
window , 1 , x s i z e=s i z , y s i z e=s i z , t i t l e = ’PHASE ’ , xpos =0 , Ypos=0
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
phase=p h a s e r e f
87
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
wset , 1
t v s c l , phase
phase=complex ( c o s ( phase ) , s i n ( phase ) )
arrow , 5 0 , 0 , 5 0 , 2 ∗mm
arrow , 0 , 5 0 ,mm, 5 0
wset , 1
;DISEGNO L ’ARRAY DI SOTTOPUPILLE
For I s u b a p =0 , N subap− 1 do begin
xpos=Centr PUP ( 0 , I s u b a p )
ypos=Centr PUP ( 1 , I s u b a p )
y1=Ypos−e x a g s i d e
y2=Ypos−e x a g s i d e /2
y3=Ypos+e x a g s i d e /2
y4=Ypos+e x a g s i d e
x1=Xpos−e x a g a p t
x2=Xpos
x3=Xpos+e x a g a p t
p l o t s , [ X2 , X3 ] , [ Y1 , Y2 ] , / d e v i c e , c o l o r =255 ; [ 0 , 0 , 0 ]
p l o t s , [ X3 , X3 ] , [ Y2 , Y3 ] , / d e v i c e , c o l o r =255 ; [ 0 , 0 , 0 ]
p l o t s , [ X3 , X2 ] , [ Y3 , Y4 ] , / d e v i c e , c o l o r =255 ; [ 0 , 0 , 0 ]
p l o t s , [ X2 , X1 ] , [ Y4 , Y3 ] , / d e v i c e , c o l o r =255 ; [ 0 , 0 , 0 ]
p l o t s , [ X1 , X1 ] , [ Y3 , Y2 ] , / d e v i c e , c o l o r =255 ; [ 0 , 0 , 0 ]
p l o t s , [ X1 , X2 ] , [ Y2 , Y1 ] , / d e v i c e , c o l o r =255 ; [ 0 , 0 , 0 ]
p l o t s , [ Xpos , Xpos ] , [ Ypos+5 ,Ypos − 5 ] , / d e v i c e , c o l o r =255 ; [ 0 , 0 , 0 ]
p l o t s , [ Xpos −5 ,Xpos + 5 ] , [ Ypos , Ypos ] , / d e v i c e , c o l o r =255 ; [ 0 , 0 , 0 ]
endfor
i d =0
; N subap− 1
For I s u b a p =0 , N subap −1 do begin
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
xpos=Centr POS ( 0 , i s u b a p )− f i x (xm / 2 . )
ypos=Centr POS ( 1 , i s u b a p )− f i x (ym / 2 . )
imm=REFERENCE( xpos : xpos+xm−1 , ypos : ypos+ym−1)
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
xpos=Centr PUP ( 0 , I s u b a p )
ypos=Centr PUP ( 1 , I s u b a p )
subpup= f l t a r r ( e x a g s i d e ∗ 2 , e x a g s i d e ∗ 2 )
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
MASK=Exag mask ∗SUPPORT( xpos−e x a g s i d e : xpos+e x a g s i d e
;
;
subpup=MASK∗PHASE( xpos−e x a g s i d e
window , 3
t v s c l , subpup
: xpos+e x a g s i d e
, ypos−e x a g s i d e : ypos+e x a g s i d e )
, ypos−e x a g s i d e : ypos+e x a g s i d e )
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
F IMM= f f t (imm, −1)
r= f f t (MASK, −1)
c r o s s= f f t ( c o n j ( r ) ∗ r , 1 )
c r o s s=c r o s s / c r o s s ( 0 , 0 )
88
P S F r ef=f l o a t ( f f t ( c r o s s , 1 ) )
P S F r ef=s h i f t ( PSF ref , e x a g s i d e , e x a g s i d e )
; c o n v o l u t i o n i n f o u r i e r s pac e
r e f= f f t (F IMM∗ c r o s s , 1 )
f= f f t ( subpup , −1)
c r o s s= f f t ( c o n j ( f ) ∗ f , 1 )
c r o s s=c r o s s / c r o s s ( 0 , 0 )
PSF=f l o a t ( f f t ( c r o s s , 1 ) )
PSF= s h i f t (PSF , e x a g s i d e , e x a g s i d e )
subap= f f t (F IMM∗ c r o s s , 1 )
; r e s u l t Ã¨ l a p s f t e o r i c a f i t t a t a con una g a u s s i a n a
R e s u l t = GAUSS2DFIT (PSF , A1 )
R e s u l t = GAUSS2DFIT ( PSF ref , A2 )
s h i f t x ( I s u b a p )=(A1(4) −A2 ( 4 ) )
s h i f t y ( I s u b a p )=(A1(5) −A2 ( 5 ) )
; print , s h i f t x ( I subap ) , s h i f t y ( I subap )
D v ect ( i d )= s h i f t x ( i s u b a p )
D v ect ( i d +1)= s h i f t y ( i s u b a p )
i d=i d +2
i f abs ( s h i f t x ( i s u b a p ) ) gt 0 . 0 or abs ( s h i f t y ( i s u b a p ) ) gt 0 . 0 then begin
t a i l =[ xpos , ypos ]
h ea d =[ xpos +50∗ s h i f t x ( I s u b a p ) , ypos +50∗ s h i f t y ( I s u b a p ) ]
wset , 1
ARROW, t a i l ( 0 ) , t a i l ( 1 ) , h ea d ( 0 ) , h ea d ( 1 ) , c o l o r =255 , / s o l i d
endif
endfor
T o t v e c t (∗ , 0 )= D v ect
Z v e c t= SVSOL(U, W, V, D v ect )
Z vect (1)=0.
TOT Z VECT(∗ , 0 )= Z v e c t / ( 2 ∗ ! p i )
z v e c t=z v e c t / ( 2 ∗ ! p i )
w d el ete , 1
d v e c t=d v e c t ∗0
z v e c t=z v e c t ∗0
mask=mask ∗0
exag mask=exag mask ∗0
s i z =500
e r r o r ( image )= t o t a l ( ( p h a s e r e f ) ˆ 2 . ) / t o t a l ( s u p p o r t )
print , ’max ZERNIKE ’
print , image
89
, h t h i c k = 2 . 0 , t h i c k =2.
; endfor
endfor
window , 4
; x=f i n d g e n ( 250)
plot , abs ( e r r o r )
; save , f i l e n a m e =’/home/ s o l a r −x2 / Desktop / e r r o r . sav ’ , e r r o r
stop
end
pro f i t t i n g e r r o r
; program t o compute Z e r n i k e f i t t i n g e r r o r up t o 200 p o l y n o m i a l o r d e r
; Marco S t a n g a l i n i a p r i l e 2008
; l a s t mod 29/04/08
; n comp Ã¨ i l numero d i componenti max da u t i l i z z a r e pe r l a f a s e d i r e f e r e n c e e l ’ andamento d
; n s t a t Ã¨ i l numero d i immagini su c u i e s e g u i r e l ’ a l g o r i t m o , de v e e s s e r e a l t o pe r a v e r e a b b a
; su q u e s t e immagini r a p p r e s e n t a n o i l numero d i r e a l i z z a z i o n i d e l l ’ e n s e m b l e
; come s t i m a d e l l ’ e r r o r e s u l f i t v i e n e p r e s o l ’ i n t e g r a l e d e l l o s c a r t o q u a d r a t i c o
; d e l l e due f a s i n o r m a l i z z a t o s u l l a s u p e r f i c i e d e l l a p u p i l l a
DEVICE, DECOMPOSED = 0
loadct ,13
n comp=200
n s t a t =2
pc=0
n s t a t=n s t a t −1
FIT ERR =d b l a r r ( 2 2 6 , n comp )
; number o f components f o r
; numero immagini e n s e m b l e
n comp=n comp−1
z i= f l t a r r ( 2 0 0 )
k i= f l t a r r ( 2 0 0 )
p h a s e0= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
p h a s e p r o z e= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
phase= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
p h a s e p r o= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
r e c z e= f l t a r r ( 2 0 0 )
r e c o v e r y= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
x=f i n d g e n ( 2 0 0 )
f i t e r r m e a n=d b l a r r ( 2 0 0 )
r e s t o r e , ’ / media/ u s b d i s k / z e r n i k e 2 0 0 G f l o a t . s a v ’ , / v e r b o s e
r e s t o r e , ’ / media/ u s b d i s k / t o t z v e c t 2 0 0 . s a v ’ , / v e r b o s e
; r e s t o r e , ’ / media / u s b d i s k / t o t z v e c t r e m e a s u r e d . sav ’ , / v e r b o s e
; se a t t i v o questas opzione bisogna
; r e s t o r e , ’ / media / u s b d i s k / t o t z v e c t p u p i l . sav ’ , / v e r b o s e
norm=1./ t o t a l ( z e r ( ∗ , ∗ , 0 ) )
t o t z v e c t (∗ ,2:3)=0
f o r immagine=1 , n s t a t do begin
p h a s e p r o=p h a s e p r o ∗0
90
p h a s e p r o z e=p h a s e p r o z e ∗0
p h a s e0=p h a s e0 ∗0
p h a s e r e f=p h a s e0 ∗0
; window , 1 , t i t l e =’ r e f ’
print , immagine
;
; REFERENCE PHASE ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
f o r z =2 ,199 do begin
p h a s e r e f=p h a s e r e f+t o t z v e c t ( z , immagine ) ∗ZER( ∗ , ∗ , z )
endfor
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
t v s c l , congrid ( phase ref ,300 ,300)
p h a s e0=p h a s e r e f
; ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ PROIEZIONE WF ON ZE ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
f o r r =2 , n comp do begin
r e c z e ( r )= t o t a l ( ( ( ZER( ∗ , ∗ , r ) ) ∗ p h a s e0 )/ t o t a l ( z e r ( ∗ , ∗ , 0 ) ) )
p h a s e0=phase0−r e c z e ( r )∗ZER( ∗ , ∗ , r )
; p r i n t , t o t z v e c t ( r , immagine ) , r e c z e ( r )
endfor
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
t=0
; window , 2
f o r z2 =2 , n comp do begin
p h a s e p r o z e=p h a s e p r o z e+r e c z e ( z2 ) ∗ZER( ∗ , ∗ , z2 )
; t=t +( ab s ( t o t z v e c t ( z2 , immagine)− r e c z e ( z2 ) ) ˆ 2 . )
f i t e r r ( immagine , z2 )=norm∗ t o t a l ( ( p h a s e r e f −p h a s e p r o z e ) ˆ 2 . )
; p r i n t , norm ∗ ( t o t a l ( z e r ( ∗ , ∗ , 0 ) ∗ p h a s e r e f −z e r ( ∗ , ∗ , 0 ) ∗ p h a s e p r o z e ) ˆ 2 )
endfor
endfor
f o r e =2 , n comp do begin
f i t e r r m e a n ( e)=mean ( f i t e r r ( ∗ , e ) )
endfor
loadct , 0
window , 2 , t i t l e = ’ phase e r r o r ’
plot , x , f i t e r r ( 1 , ∗ ) , psym=2 , x t i t l e= ’ Z e r n i k e mode number ’ , y t i t l e= ’ phase r e c o n s t r u c t i o n
e r r o r [ waves ] ’
print , z2
r e s t o r e , ’ / media / u s b d i s k / f i t t i n g e r r o r z e m e a s u r e d . s a v ’ , / v e r b o s e
loadct , 3
o p l o t , f i t e r r m e a s u r e d ( 0 , ∗ ) , c o l o r =100 , psym=5
; save , f i l e n a m e =’/media / u s b d i s k / p h a s e r e f . sav ’ , p h a s e r e f
; save , f i l e n a m e =’/media / u s b d i s k / f i t t i n g e r r o r z e m e a s u r e d . sav ’ , f i t e r r
stop
end
pro f i t t i n g e r r o r k l
; program t o compute Karhunen−Loeve f i t t i n g e r r o r up t o 200 p o l y n o m i a l o r d e r
; Marco S t a n g a l i n i a p r i l e 2008
; l a s t mod 30/04/08
; n comp Ã¨ i l numero d i componenti max da u t i l i z z a r e pe r l a f a s e d i r e f e r e n c e e l ’ andamento d e l
91
fittin
; n s t a t Ã¨ i l numero d i immagini su c u i e s e g u i r e l ’ a l g o r i t m o , de v e e s s e r e a l t o pe r a v e r e a b b a
; su q u e s t e immagini r a p p r e s e n t a n o i l numero d i r e a l i z z a z i o n i d e l l ’ e n s e m b l e
; come s t i m a d e l l ’ e r r o r e s u l f i t v i e n e p r e s o l ’ i n t e g r a l e d e l l o s c a r t o q u a d r a t i c o
; d e l l e due f a s i n o r m a l i z z a t o s u l l a s u p e r f i c i e d e l l a p u p i l l a
DEVICE, DECOMPOSED = 0
loadct ,13
n comp=200
n s t a t =2
pc=0
n s t a t=n s t a t −1
FIT ERR = f l t a r r ( 2 2 6 , n comp )
n comp=n comp−1
z i= f l t a r r ( 2 0 0 )
k i= f l t a r r ( 2 0 0 )
p h a s e0= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
p h a s e p r o k l= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
phase= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
p h a s e p r o= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
r e c k l= f l t a r r ( 2 0 0 )
c j= f l t a r r ( 2 0 0 )
r e c o v e r y= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
x=f i n d g e n ( 2 0 0 )
f i t e r r m e a n k l=d b l a r r ( 2 0 0 )
; number o f components f o r
; numero immagini e n s e m b l e
r e s t o r e , ’ / media/ u s b d i s k / K L 2 0 0 f l o a t . s a v ’ , / v e r b o s e
r e s t o r e , ’ / media/ u s b d i s k / t o t z v e c t 2 0 0 . s a v ’ , / v e r b o s e
r e s t o r e , ’ / media/ u s b d i s k / K L 2 0 0 e i g e n v e c t . s a v ’ , / v e r b o s e
r e s t o r e , ’ / media/ u s b d i s k / z e r n i k e 2 0 0 G f l o a t . s a v ’ , / v e r b o s e
norm=1./ t o t a l ( z e r ( ∗ , ∗ , 0 ) )
; t o t z v e c t ( ∗ , 2: 3) = 0
FOR immagine=1 , n s t a t DO BEGIN
p h a s e p r o=p h a s e p r o ∗0
p h a s e p r o k l=p h a s e p r o k l ∗0
p h a s e0=p h a s e0 ∗0
p h a s e r e f=p h a s e0 ∗0
; window , 1 , t i t l e =’ r e f ’
print , immagine
; z i 0=t o t z v e c t ( 4 : n comp , immagine )
; ;A=( z i 0 # ( t r a n s p o s e ( z i 0 ) ) )
;B =e i g e n v e c t # z i 0
;
; c j [ 0 : 3 ] = t o t z v e c t ( 0 : 3 , immagine )
; c j [ 4 : n comp ]= z i 0
; r e s t o r e , ’ / media / u s b d i s k / p h a s e r e f . sav ’
; REFERENCE PHASE ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
f o r z =2 ,199 do begin
p h a s e r e f=p h a s e r e f+t o t z v e c t ( z , immagine ) ∗ZER( ∗ , ∗ , z )
endfor
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
;
; REFERENCE PHASE ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
;
f o r z =2 ,199 do b e g i n
92
;
;
;
;;
p h a s e r e f=p h a s e r e f+c j ( z ) ∗ karhunen ( ∗ , ∗ , z )
endfor
;∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
t v s c l , congrid ( phase ref ,300 ,300)
p h a s e0=p h a s e r e f
; ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ PROIEZIONE WF ON ZE ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
;
f o r r =n comp ,2 , −1 do b e g i n
f o r r =2 , n comp do begin
i f r l e 3 then begin
r e c k l ( r )= t o t a l ( ( ( Zer ( ∗ , ∗ , r ) ) ∗ p h a s e0 )/ t o t a l ( z e r ( ∗ , ∗ , 0 ) ) )
p h a s e0=phase0−r e c k l ( r )∗ z e r ( ∗ , ∗ , r )
endif e l s e begin
r e c k l ( r )= t o t a l ( ( ( karhunen ( ∗ , ∗ , r ) ) ∗ p h a s e0 )/ t o t a l ( z e r ( ∗ , ∗ , 0 ) ) )
p h a s e0=phase0−r e c k l ( r )∗ karhunen ( ∗ , ∗ , r )
; p r i n t , t o t z v e c t ( r , immagine ) , r e c z e ( r )
endelse
endfor
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
f o r t =0 ,3 do begin
p h a s e p r o k l=p h a s e p r o k l+r e c k l ( t )∗ z e r ( ∗ , ∗ , t )
f i t e r r ( immagine , t )=norm∗ t o t a l ( ( p h a s e r e f −p h a s e p r o k l ) ˆ 2 . )
stop
endfor
f o r z2=n comp ,4 , −1 do begin
; f o r z2 =2,n comp do b e g i n
;
;
;
;
i f z2 l e 3 t he n b e g i n
p h a s e p r o k l=p h a s e p r o k l+r e c k l ( z2 ) ∗ z e r ( ∗ , ∗ , z2 )
f i t e r r ( immagine , z2)=norm∗ t o t a l ( ( p h a s e r e f −p h a s e p r o k l ) ˆ 2 . )
endif else begin
p h a s e p r o k l=p h a s e p r o k l+r e c k l ( z2 )∗ karhunen ( ∗ , ∗ , z2 )
f i t e r r ( immagine , z2 )=norm∗ t o t a l ( ( p h a s e r e f −p h a s e p r o k l ) ˆ 2 . )
; endelse
; stop
; p r i n t , f i t e r r ( immagine , z2 )
stop
endfor
ENDFOR
; f o r e =2,n comp do b e g i n
; f i t e r r m e a n k l ( e)=mean ( f i t e r r ( ∗ , e ) )
; endfor
loadct , 0
window , 0 , t i t l e = ’ phase e r r o r ’
f i t e r r m e a n k l (∗ )= ( f i t e r r ( 1 , ∗ ) )
f i t e r r m e a n k l (4:199)= r ever s e ( f i t e r r m e a n k l ( 4 : 1 9 9 ) )
plot , x , f i t e r r m e a n k l , psym=2 , x t i t l e= ’KL mode number ’ , y t i t l e= ’ phase r e c o n s t r u c t i o n
93
e r r o r [ waves ] ’
; save , f i l e n a m e =’/media / u s b d i s k / f i t t i n g e r r o r k l . sav ’ , f i t e r r m e a n k l
stop
end
pro f i t t i n g e r r o r M I
r e s t o r e , ’ / media/ d i s k / karhunen MI . s a v ’
; r e s t o r e , ’ / media / d i s k / phaseMI . sav ’
; r e s t o r e , ’ / media / d i s k / K L 2 0 0 f l o a t . sav ’
r e s t o r e , ’ / media/ d i s k / t o t z v e c t 2 0 0 . s a v ’ , / v e r b o s e
r e s t o r e , ’ / media/ d i s k / z e r n i k e 2 0 0 G f l o a t . s a v ’ , / v e r b o s e
norm=1./ t o t a l ( z e r ( ∗ , ∗ , 0 ) )
p h a s e0= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
p h a s e p r o k l= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
phase= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
p h a s e p r o= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
; k l=karhunen
; k l (4:199)= r e v e r s e ( karhunen ( 4 : 1 9 9 ) )
r e c k l= f l t a r r ( 2 0 0 )
r e c= f l t a r r ( 2 0 0 )
f i t e r r=f l t a r r (200)
immagine=1
p h a s e p r o=p h a s e p r o ∗0
p h a s e p r o k l=p h a s e p r o k l ∗0
p h a s e0=p h a s e0 ∗0
p h a s e r e f=p h a s e0 ∗0
f o r z =2 ,199 do begin
p h a s e r e f=p h a s e r e f+t o t z v e c t ( z , immagine ) ∗ZER( ∗ , ∗ , z )
endfor
p h a s e0=p h a s e r e f
f o r r =2 ,199 do begin
i f r l e 3 then begin
r e c k l ( r )= t o t a l ( ( ( Zer ( ∗ , ∗ , r ) ) ∗ p h a s e0 )/ t o t a l ( z e r ( ∗ , ∗ , 0 ) ) )
p h a s e0=phase0−r e c k l ( r )∗ z e r ( ∗ , ∗ , r )
endif e l s e begin
r e c k l ( r )= t o t a l ( ( ( k l ( ∗ , ∗ , r ) ) ∗ p h a s e0 )/ t o t a l ( z e r ( ∗ , ∗ , 0 ) ) )
p h a s e0=phase0−r e c k l ( r )∗ k l ( ∗ , ∗ , r )
endelse
endfor
p h a s e r e f=p h a s e r e f ∗ 0 .
z=0
f o r z =0 ,3 do begin
p h a s e r e f=p h a s e r e f+r e c k l ( z )∗ z e r ( ∗ , ∗ , z )
endfor
94
f o r z =4 ,199 do begin
p h a s e r e f=p h a s e r e f+r e c k l ( z )∗ k l ( ∗ , ∗ , z )
endfor
f o r t =0 ,3 do begin
p h a s e p r o k l=p h a s e p r o k l+r e c k l ( t )∗ z e r ( ∗ , ∗ , t )
f i t e r r ( t )=norm∗ t o t a l ( ( p h a s e r e f −p h a s e p r o k l ) ˆ 2 . )
endfor
f o r z2 =199 ,4 , −1 do begin
; f o r z2 =4 ,199 do b e g i n
p h a s e p r o k l=p h a s e p r o k l+r e c k l ( z2 )∗ k l ( ∗ , ∗ , z2 )
f i t e r r ( z2 )=norm∗ t o t a l ( ( p h a s e r e f −p h a s e p r o k l ) ˆ 2 . )
endfor
f i t e r r (4:199)= r ever s e ( f i t e r r ( 4 : 1 9 9 ) )
;FOR I =199,0, −1 DO BEGIN
;
REC( I )=TOTAL( ( (KL( ∗ , ∗ , I ) ) ∗PHASE REF)/TOTAL(KL( ∗ , ∗ , 0 ) ) )
;
PHASE REF=PHASE REF−REC( I ) ∗KL( ∗ , ∗ , I )
;
STOP
;ENDFOR
; f o r t =0 ,199 do b e g i n
;
phas e 0=phas e 0+r e c ( t )∗ k l ( ∗ , ∗ , t )
; endfor
;
;
; f o r z =0 ,199 do b e g i n
;
p h a s e p r o k l=p h a s e p r o k l+r e c ( z ) ∗ k l ( ∗ , ∗ , z )
;
f i t e r r ( z )=norm∗ t o t a l ( ( phase0−p h a s e p r o k l ) ˆ 2 . )
; endfor
loadct , 0
plot , f i t e r r
stop
end
pro f i t t i n g e r r o r M I
r e s t o r e , ’ / media / d i s k / karhunen MI . s a v ’
; r e s t o r e , ’ / media / d i s k /phaseMI . sav ’
; r e s t o r e , ’ / media / d i s k / K L 2 0 0 f l o a t . sav ’
r e s t o r e , ’ / media / d i s k / t o t z v e c t 2 0 0 . s a v ’ , / v e r b o s e
r e s t o r e , ’ / media / d i s k / z e r n i k e 2 0 0 G f l o a t . s a v ’ , / v e r b o s e
norm=1./ t o t a l ( z e r ( ∗ , ∗ , 0 ) )
p h a s e0= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
p h a s e p r o k l= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
phase= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
p h a s e p r o= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
; k l=karhunen
; k l (4:199)= r e v e r s e ( karhunen ( 4 : 1 9 9 ) )
r e c k l= f l t a r r ( 2 0 0 )
r e c= f l t a r r ( 2 0 0 )
95
f i t e r r=f l t a r r (200)
immagine=1
p h a s e p r o=p h a s e p r o ∗0
p h a s e p r o k l=p h a s e p r o k l ∗0
p h a s e0=p h a s e0 ∗0
p h a s e r e f=p h a s e0 ∗0
f o r z =2 ,199 do begin
p h a s e r e f=p h a s e r e f+t o t z v e c t ( z , immagine ) ∗ZER( ∗ , ∗ , z )
endfor
p h a s e0=p h a s e r e f
f o r r =2 ,199 do begin
i f r l e 3 then begin
r e c k l ( r )= t o t a l ( ( ( Zer ( ∗ , ∗ , r ) ) ∗ p h a s e0 )/ t o t a l ( z e r ( ∗ , ∗ , 0 ) ) )
p h a s e0=phase0−r e c k l ( r )∗ z e r ( ∗ , ∗ , r )
endif e l s e begin
r e c k l ( r )= t o t a l ( ( ( k l ( ∗ , ∗ , r ) ) ∗ p h a s e0 )/ t o t a l ( z e r ( ∗ , ∗ , 0 ) ) )
p h a s e0=phase0−r e c k l ( r )∗ k l ( ∗ , ∗ , r )
endelse
endfor
p h a s e r e f=p h a s e r e f ∗ 0 .
z=0
f o r z =0 ,3 do begin
p h a s e r e f=p h a s e r e f+r e c k l ( z )∗ z e r ( ∗ , ∗ , z )
endfor
f o r z =4 ,199 do begin
p h a s e r e f=p h a s e r e f+r e c k l ( z )∗ k l ( ∗ , ∗ , z )
endfor
f o r t =0 ,3 do begin
p h a s e p r o k l=p h a s e p r o k l+r e c k l ( t )∗ z e r ( ∗ , ∗ , t )
f i t e r r ( t )=norm∗ t o t a l ( ( p h a s e r e f −p h a s e p r o k l ) ˆ 2 . )
endfor
f o r z2 =199 ,4 , −1 do begin
; f o r z2 =4 ,199 do b e g i n
p h a s e p r o k l=p h a s e p r o k l+r e c k l ( z2 )∗ k l ( ∗ , ∗ , z2 )
f i t e r r ( z2 )=norm∗ t o t a l ( ( p h a s e r e f −p h a s e p r o k l ) ˆ 2 . )
endfor
f i t e r r (4:199)= r ever s e ( f i t e r r ( 4 : 1 9 9 ) )
;FOR I =199,0, −1 DO BEGIN
;
REC( I )=TOTAL( ( (KL( ∗ , ∗ , I ) ) ∗PHASE REF)/TOTAL(KL( ∗ , ∗ , 0 ) ) )
;
PHASE REF=PHASE REF−REC( I ) ∗KL( ∗ , ∗ , I )
;
STOP
;ENDFOR
96
; f o r t =0 ,199 do b e g i n
;
phas e 0=phas e 0+r e c ( t )∗ k l ( ∗ , ∗ , t )
; endfor
;
;
; f o r z =0 ,199 do b e g i n
;
p h a s e p r o k l=p h a s e p r o k l+r e c ( z ) ∗ k l ( ∗ , ∗ , z )
;
f i t e r r ( z )=norm∗ t o t a l ( ( phase0−p h a s e p r o k l ) ˆ 2 . )
; endfor
loadct , 0
plot , f i t e r r
stop
end
pro f i t t i n g e r r o r r e m e a s u r e d
; program t o compute Z e r n i k e f i t t i n g e r r o r up t o 200 p o l y n o m i a l o r d e r
; Marco S t a n g a l i n i a p r i l e 2008
; l a s t mod 29/04/08
; n comp Ã¨ i l numero d i componenti max da u t i l i z z a r e pe r l a f a s e d i r e f e r e n c e e l ’ andamento d e l f i t t i n
; n s t a t Ã¨ i l numero d i immagini su c u i e s e g u i r e l ’ a l g o r i t m o , de v e e s s e r e a l t o pe r a v e r e a b b a s t a n z a s
; su q u e s t e immagini r a p p r e s e n t a n o i l numero d i r e a l i z z a z i o n i d e l l ’ e n s e m b l e
; come s t i m a d e l l ’ e r r o r e s u l f i t v i e n e p r e s o l ’ i n t e g r a l e d e l l o s c a r t o q u a d r a t i c o
; d e l l e due f a s i n o r m a l i z z a t o s u l l a s u p e r f i c i e d e l l a p u p i l l a
DEVICE, DECOMPOSED = 0
loadct ,13
n comp=200
n s t a t =2
pc=0
n s t a t=n s t a t −1
FIT ERR =d b l a r r ( 2 2 6 , n comp )
; number o f components f o r PCA
; numero immagini e n s e m b l e
n comp=n comp−1
z i= f l t a r r ( 2 0 0 )
k i= f l t a r r ( 2 0 0 )
p h a s e0= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
p h a s e p r o z e= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
phase= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
p h a s e p r o= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
r e c z e= f l t a r r ( 2 0 0 )
r e c o v e r y= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
x=f i n d g e n ( 2 0 0 )
f i t e r r m e a n=d b l a r r ( 2 0 0 )
r e s t o r e , ’ / media / u s b d i s k / z e r n i k e 2 0 0 G f l o a t . s a v ’ , / v e r b o s e
r e s t o r e , ’ / media / u s b d i s k / t o t z v e c t r e m e a s u r e d . s a v ’ , / v e r b o s e
r e s t o r e , ’ / media / u s b d i s k / p h a s e r e f . s a v ’ , / v e r b o s e
; se a t t i v o questas opzione bisogna s o s t i t u i r e
norm=1./ t o t a l ( z e r ( ∗ , ∗ , 0 ) )
; t o t z v e c t ( ∗ , 2: 3) = 0
; p h a s e r e f=phas e 0 ∗0
p h a s e p r o=p h a s e p r o ∗0
p h a s e p r o z e=p h a s e p r o z e ∗0
97
tot
p h a s e0=p h a s e0 ∗0
; window , 1 , t i t l e =’ r e f ’
;
;
; REFERENCE PHASE ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
;
f o r z =2 ,199 do b e g i n
;
p h a s e r e f=p h a s e r e f+z v e c t ( z ) ∗ZER( ∗ , ∗ , z )
;
endfor
;
;∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
;
t v s c l , congrid ( phase ref ,300 ,300)
;
; phas e 0=p h a s e r e f
; ; ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ PROIEZIONE WF ON ZE ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
;
f o r r =2,n comp do b e g i n
;
r e c z e ( r)= t o t a l ( ( (ZER( ∗ , ∗ , r ) ) ∗ phas e 0 )/ t o t a l ( z e r ( ∗ , ∗ , 0 ) ) )
;
phas e 0=phase0−r e c z e ( r ) ∗ZER( ∗ , ∗ , r )
;
; p r i n t , t o t z v e c t ( r , immagine ) , r e c z e ( r )
;
endfor
; ;∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
t=0
; window , 2
f o r z2 =2 , n comp do begin
p h a s e p r o z e=p h a s e p r o z e+z v e c t ( z2 )∗ZER( ∗ , ∗ , z2 )
f i t e r r ( 0 , z2 )=norm∗ t o t a l ( ( p h a s e r e f −p h a s e p r o z e ) ˆ 2 . )
endfor
f o r e =2 , n comp do begin
f i t e r r m e a n ( e)=mean ( f i t e r r ( ∗ , e ) )
endfor
loadct , 0
window , 2 , t i t l e = ’ phase e r r o r ’
plot , x , f i t e r r ( 0 , ∗ ) , psym=2 , x t i t l e= ’ Z e r n i k e mode number ’ , y t i t l e= ’ phase r e c o n s t r u c t i o n
error [ w
f i t e r r m e a s u r e d= f i t e r r
print , z2
; save , f i l e n a m e =’/media / u s b d i s k / p h a s e r e f . sav ’ , p h a s e r e f
save , f i l e n a m e= ’ / media/ u s b d i s k / f i t t i n g e r r o r z e m e a s u r e d . s a v ’ , f i t e r r m e a s u r e d
stop
end
pro k a r h u n e n a n n u l a r
; This program d e f i n e s Karhunen−Loeve p o l y n o m i a l b a s i s and o p e r a t e s
; phas e p r o j e c t i o n on i t comparing t h a t w i t h Z e r n i k e r e c o n s t r u c t i o n
; Marco S t a n g a l i n i 2008
DEVICE, DECOMPOSED = 0
loadct ,13
r e s t o r e , ’ / media/ u s b d i s k / t o t z v e c t 2 0 0 . s a v ’ , / v e r b o s e
r e s t o r e , ’ / media/ u s b d i s k / z e r n i k e 2 0 0 G f l o a t . s a v ’ , r e s t o r e d o b j e c t s=ZE , / v e r b o s e
immagine=203
; wor k i n g image
; numero immagini da mediare p
n s t a t =2
n c o e f f =200
98
n s t a t=n s t a t −1
p h a s e2= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
p h a s e0= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
Karhunen= f l t a r r ( 1 0 4 4 , 1 0 4 4 , n c o e f f )
phase= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
U= f l t a r r ( n c o e f f , n c o e f f )
kk= f l t a r r ( 1 9 6 )
window , 0 , t i t l e = ’ phase ZE ’ , x s i z e =300 , y s i z e =300
z i= f l t a r r ( n c o e f f )
;LEARNING
f o r u=0 , n c o e f f −1 do begin
z i ( u)=mean ( t o t z v e c t ( u , 0 : n s t a t ) )
endfor
; S e l e z i o n o s o l o i c o e f f i c i e n t i >3 i n modo da f i s s a r e t i p − t i l t con Z e r n i k e
z i 0=z i ( 4 : n c o e f f −1)
z i 0=z i 0 −mean ( z i 0 )
A=( z i 0 # ( t r a n s p o s e ( z i 0 ) ) )
;
; Calcolo d e l l a camtrice contenente g l i a u t o v e t t o r i d e l l a covarianza
;AUTOVALORI E AUTOVETTORI
t r i r e d , A, D, E
; D c o n t a i n s d i a g o n a l , E c o n t a i n s o f f −d i a g o n a l
t r i q l , D, E , A
; D c o n t a i n s t h e e i g e n −v a l u e s , U = e i g e n v e c t o r s
e i g e n v a l = (D)
e i g e n v e c t =t r a n s p o s e (A)
;!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
print , s i z e (A)
i n d e x = s o r t (D)
e i g e n v e c t =t r a n s p o s e (A)
; Order by i n c r e a s i n g e i g e n v a l u e
e i g e n v e c t=e i g e n v e c t [ index , ∗ ]
;
;
;
e i g e n v a l = HQR(ELMHES(A) , /DOUBLE)
e i g e n v e c t = EIGENVEC(A, e i g e n v a l )
; CALCOLO COEFFICIENTI ESPANSIONE KARHUNEN−LOEVE
; The c o e f f i c i e n t s b j o f t h e K−L e x pan s i on ar e t h e components o f t h e column v e c t o r B=U∗ZE
B =r e a l p a r t ( e i g e n v e c t ) # z i 0
stop
; ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ NORMALIZZZIONE ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
; f o r c o l =0,16 do b e g i n
; e i g e n v e c t ( c o l ,∗)= e i g e n v e c t ( c o l , ∗ ) / t o t a l ( e i g e n v e c t ( c o l , ∗ ) )
; endfor
; wavefront z e r n i k e
f o r z =1 , n c o e f f −1 do begin
phase=phase+z i ( z ) ∗ z e r ( ∗ , ∗ , z )
endfor
wset , 0
t v s c l , c o n g r i d ( phase , 3 0 0 , 3 0 0 )
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
; c a l c o l o l ’ i n v e r s a pe r e s p r i m e r e i c o e f f d i z e r n i k e
U=t r a n s p o s e (A)
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
99
; f o r z o i pr i m i t r e p o l i n o m i ad e s s e r e u g u a l i a q u e l l i d i z e r n i k e
f o r k =0 ,3 do begin
karhunen ( ∗ , ∗ , k)= z e r ( ∗ , ∗ , k )
endfor
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
; ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ BASE KARHUNEN−LOEVE ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
; S c r i v o i KL i n f u n z i o n e d e g l i Z e r n i k e con una sommatoria
f o r t =4 , n c o e f f −1 do begin
Karhunen ( ∗ , ∗ , t )=0.
f o r r =4 , n c o e f f −1 do Karhunen ( ∗ , ∗ , t )=U( t −4 , r −4)∗ZER( ∗ , ∗ , r )+
endfor
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
; w a v e f r o n t Karhunen
f o r z =1 , n c o e f f −1 do begin
i f z l e 3 then begin
p h a s e2=p h a s e2+z i ( z ) ∗ZER( ∗ , ∗ , z )
endif e l s e begin
p h a s e2=p h a s e2+B( z −4)∗ Karhunen ( ∗ , ∗ , z )
endelse
endfor
window , 2 , x s i z e =300 , y s i z e =300
t v s c l , c o n g r i d ( phase2 , 3 0 0 , 3 0 0 )
print , ’max r e s ’
print , max( phase−p h a s e2 )
print , ’ s t a t i s t i c s ’
r e s u l t = MOMENT( phase−p h a s e2 )
PRINT, ’ Mean : ’ , r e s u l t [ 0 ] & PRINT, ’ V a r i a n ce : ’ , r e s u l t [ 1 ] & $
PRINT, ’ Skewness : ’ , r e s u l t [ 2 ] & PRINT, ’ K u r t o s i s : ’ , r e s u l t [ 3 ]
; save , f i l e n a m e =’/media / u s b d i s k / K L 200 e i g e n v e c t . sav ’ , e i g e n v e c t
; save , f i l e n a m e =’/media / u s b d i s k / K L 2 0 0 f l o a t . sav ’ , karhunen
stop
end
pro LP
; t h i s program computes ARMA( p , q ) p r o c e s s par am e t e r s and
; t ry to f o r e c a s t values of Zernike c o e f f i c i e n t s
print , s y s t i m e ( )
DEVICE, DECOMPOSED = 0
100
d a t i r e a l i =0
; 0=OFF 1=ON
i f d a t i r e a l i eq 0 then begin
r e s t o r e , ’ / media / d i s k −3/ t o t z v e c t t s 2 . s a v ’ , / v e r b o s e
t o t z v e c t t s=t o t z v e c t t s ∗ 1 0 .
t o t z v e c t t s (∗ , 1 7 )= t o t z v e c t t s ( ∗ , 1 6 )
endif e l s e begin
r e s t o r e , ’ / media / d i s k −3/ t o t z v e c t 2 0 0 . s a v ’ , / v e r b o s e
t o t z v e c t t s=t o t z v e c t
endelse
; r e s t o r e , ’ / media / d i s k / t o t z v e c t t s . sav ’ , / v e r b o s e
; t o t z v e c t t s (∗ ,17)= t o t z v e c t t s ( ∗ , 1 6 )
; t o t z v e c t t s =t o t z v e c t
; t o t z v e c t t s =t o t z v e c t t s ∗ 1 0 .
max ze=35
x= f l t a r r ( 2 0 0 )
lagmax=100
l e a r n i n g =100
e s t i m a t e= f l t a r r ( 2 0 0 )
p=100
q=0
time max =100
array AR= f l t a r r ( p )
; array MA= f l t a r r ( q )
z= f l t a r r ( 1 0 0 )
e s t=e s t i m a t e ∗0
s t i m e= f l t a r r ( 5 0 0 )
ze =16
; f o r qq =0,50 do b e g i n
; f o r z e =2 ,199 do b e g i n
;
f o r t=0+qq , l e a r n i n g+qq do b e g i n
f o r t =0 , l e a r n i n g do begin
x ( t )= t o t z v e c t t s ( ze , t )
endfor
; i f q g t 0 t he n b e g i n
; x ( l e a r n i n g+qq+1)= e s t i m a t e ( z e )
; endif
t=0
ACF=IMSL AUTOCORRELATION( x , lagmax )
r o=a c f / a c f ( 0 )
mu=mean ( x )
; window , 0
; p l o t , ac f , x t i t l e =’Time Lag (1 s t e p = 10 ms ) ’ , y t i t l e =’ACF’
; SET PLOT, ’PS ’
; DEVICE, FILE=’/home/ s o l a r −x2 / Desktop /ACF sim . ps ’ , /COLOR, BITS=8
; p l o t , ac f , x t i t l e =’Time Lag (1 s t e p = 10 ms ) ’ , y t i t l e =’ACF’
; device ,/ c los e
; set plot , ’x ’
101
f o r e c a s t=i m s l a r m a ( x , p , q , a u to co v=a c f 1 ) ; , f o r e c a s t=f o r e , n p r e d i c t =12, B ac k war d O r i g i n = 3)
; , INIT EST AR=array AR , INIT EST MA=array MA )
ACF est=a c f 1 / a c f 1 ( 0 )
e s t i m a t e ( ze )= t o t z v e c t t s ( ze , l e a r n i n g ) ∗ f o r e c a s t (1)+ t o t z v e c t t s ( ze , l e a r n i n g −1)∗ f o r e c a s t (2)+
; r e v=r e v e r s e ( t o t z v e c t t s ( ∗ , 0 : l e a r n i n g ) , 2 )
print , ’ p r e v i s i o n e ’ , t o t z v e c t t s ( ze , l e a r n i n g ) ∗ f o r e c a s t (1)+ t o t z v e c t t s ( ze , l e a r n i n g −1)∗ f o r e c a
print , ’ r e a l e ’ , t o t z v e c t t s ( ze , l e a r n i n g +1)
; coe=
; e s t=f o r e c a s t ( 1 : p )∗ r e v ( ze , ∗ )
f o r u=1 ,p do begin
e s t ( ze )= e s t ( ze )+ f o r e c a s t ( u ) ∗ t o t z v e c t t s ( ze , l e a r n i n g −u+1)
endfor
plot , a c f
o p l o t , a c f e s t , psym=2
stop
; SET PLOT, ’PS ’
; DEVICE, FILE=’/home/ s o l a r −x2 / Desktop /ACF con . ps ’ , /COLOR, BITS=8
; ! P .MULTI = [ 0 , 2 , 1 ]
; p l o t , ac f , x t i t l e =’Time Lag (1 s t e p = 10 ms ) ’ , y t i t l e =’ACF’
; p l o t , a c f e s t , x t i t l e =’Time Lag (1 s t e p = 10 ms ) ’ , y t i t l e =’ACF’ , x r an g e = [ 0 , 100]
; device ,/ c los e
; set plot , ’ x ’
; print , ’ best estimation
’ , e st ( ze )
; endfor
;
; SET PLOT, ’PS ’
; DEVICE, FILE=’/home/ s o l a r −x2 / Desktop / f o r e c a s t r e a l . ps ’ , /COLOR, BITS=8
; ! P .MULTI = [ 0 , 1 , 3 ]
; ; window , 0 , x s i z e =1100 , y s i z e =600, t i t l e =’ Z e r n i k e f o r e c a s t i n g ’
; p l o t , t o t z v e c t t s ( 0 : 3 5 , l e a r n i n g +1) , x t i t l e =’ Z e r n i k e c o e f f i c i e n t ’ , y t i t l e =’Measured v a l u e [ wav
; o p l o t , t o t z v e c t t s ( 0 : 3 5 , l e a r n i n g +1) ,psym=5
;;
; ; ; loadct ,3
; p l o t , e s t i m a t e ( 0 : 3 5 ) , x t i t l e =’ Z e r n i k e c o e f f i c i e n t ’ , y t i t l e =’1 ms F o r e c a s t i n g [ waves ] ’
; o p l o t , e s t i m a t e ( 0 : 3 5 ) , psym=5
;;
;;
; ; p l o t , ( ( t o t z v e c t t s ( 0 : 3 5 , l e a r n i n g +10)− e s t i m a t e ) ˆ ( 2 . ) ) ˆ ( 0 . 5 ) , y r an g e = [ 0 , 0 . 7 ] , x t i t l e =’ Z e r n i k e
; ; ; loadct ,3
; p l o t , ( ( t o t z v e c t t s ( 0 : 3 5 , l e a r n i n g+1)− e s t ) ˆ ( 2 . ) ) ˆ ( 0 . 5 ) , psym=2, y r an g e = [ 0 , 0 . 7 ] , x t i t l e =’ Z e r n i k e
; device ,/ c los e
; set plot , ’ x ’
; s t i m e ( l e a r n i n g+qq)= e s t i m a t e ( z e )
! p . m u l t i=0
SET PLOT, ’ PS ’
DEVICE, FILE= ’ / media / d i s k / t e s i M a r c o / immagini / errorACF . ps ’ , / c o l o r , b i t s =8
; p l o t , ACF est , x t i t l e =’Lag ’ , y t i t l e =’ACF’ , psym=4, x r an g e = [ 0 , 100]
plot , ( ( a c f −a c f e s t ) ˆ 2 . ) ˆ 0 . 5 , y r a n g e = [ 0 , 0 . 0 0 0 0 0 0 0 8 ] , psym=2 , y t i t l e= ’ E s ti m a ted ACF a b s o l u t e r e s i d
device , / c l o s e
set plot , ’x ’
stop
print , s y s t i m e ( )
s t i m e ( 0 : l e a r n i n g )= t o t z v e c t t s ( ze , 0 : l e a r n i n g )
102
window , 2 , x s i z e =500 , y s i z e =500
; p l o t , s t i m e , x r an g e = [ 0 : 2 0 0 ]
plot , t o t z v e c t t s ( ze , ∗ ) , x r a n g e = [ 1 0 0 , 1 1 0 ]
o p l o t , s ti m e , psym=5
loadct , 3
o p l o t , s ti m e , psym=5 , c o l o r =200
stop
end
pro MTBS
;MTBS ( M u l t i l a y e r TurBulence S i m u l a t o r ) i s a s i m u l a t i o n code d e v e l o p e d f o r
; a t m o s p h e r i c t u r b o l e n c e and A dapt i v e O p t i c s d e s i g n s t u d y .
; I t can perform m u l t i l a y e r s i m u l a t i o n w i t h t h e p o s s i b i l i t y t o c hoos e
; h e i g h t and wind s pe e d . Phase s c r e e n ar e c a l c u l a t e d u s i n g PSG r o u t i n e .
; o t h e r u s e f u l par am e t e r s t o s e t u p ar e : d e t e c t o r p i x e l s c a l e and t i m e s t e p f o r a c q u i s i t i o n
;
;
;
;
;
;
;
;
;
;
;
;
This program i s f r e e s o f t w a r e : you can r e d i s t r i b u t e i t and/ or modify
i t under t h e terms o f t h e GNU General P u b l i c L i c e n s e as p u b l i s h e d by
t h e Free S o f t w a r e Foundation , e i t h e r v e r s i o n 3 o f t h e L i c e ns e , or
( a t your o p t i o n ) any l a t e r v e r s i o n .
This program i s d i s t r i b u t e d i n t h e hope t h a t i t w i l l be u s e f u l ,
b u t WITHOUT ANY WARRANTY; w i t h o u t even t h e i m p l i e d war r an t y o f
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See t h e
GNU General P u b l i c L i c e n s e f o r more d e t a i l s .
You s h o u l d have r e c e i v e d a copy o f t h e GNU General P u b l i c L i c e n s e
a l o n g w i t h t h i s program .
I f not , s e e <h t t p : / /www . gnu . or g / l i c e n s e s />.
; Author : Marco S t a n g a l i n i ( f i r s t v e r s i o n 5 j u n e 2008)
t i m e s t e p =10.
p i x s c a l e =0.2
L 1 s p eed =2.
L 2 s p eed =30.
h1 =10.
h2 =300.
l 1 =4.
l 2 =4.
; t i m e s t e p i n ms
p i x e l scale arcsec/ pix
s pe e d Layer 1 i n m/ s e c
s pe e d Layer 2 i n m/ s e c
h e i g t h in meters l a y e r 1
; h e i g t h in meters l a y e r 2
;
;
;
;
DEVICE, DECOMPOSED = 0
loadct ,13
r e s t o r e , ’ / media / d i s k / s u p p o r t . s a v ’
r e s t o r e , ’ / media / d i s k / s cr een G . s a v ’ , / v e r b o s e
s cr een G=smooth ( s c r e e n , 2 0 )
s c r e e n=s c r e e n ∗0
r e s t o r e , ’ / media / d i s k / s c r e e n L . s a v ’ , / v e r b o s e
s c r e e n L=smooth ( s c r e e n , 2 0 )
Npix =2088
rad =( p i x s c a l e ∗ ! p i ) / ( 3 6 0 0 . ∗ 1 8 0 . )
; d e t e c t o r p i x e l s c a l e ov e r t h e l a y e r
s c a l e L 1=rad ∗ h1
s c a l e L 2=rad ∗ h2
103
L 1 s p eed m s=L 1 s p eed ∗ ( 1 0 ˆ ( − 6 . ) )
L 2 s p eed m s=L 2 s p eed ∗ ( 1 0 ˆ ( − 6 . ) )
; l a y e r d i s p l a c e m e n t i n t h e t i m e s t e p ab ov e
L 1 s t e p=L 1 s p eed m s ∗ t i m e s t e p
L 2 s t e p=L 2 s p eed m s ∗ t i m e s t e p
; ; angular displacement in time s t e p
; D t he t a L 1=L 1 s t e p /h1
; D t he t a L 2=L 2 s t e p /h2
;m pe r p i x e l
p i x l 1=l 1 / Npix
p i x l 2=l 2 / Npix
; e q u i v a l e n t p i x e l displacement in t he time s t e p
Nstep L1=L 1 s t e p / p i x l 1
Nstep L2=L 2 s t e p / p i x l 2
; d e t e c t o r p i x e l s c a l e i n p i x e l ov e r t h e l a y e r
p p i x l 1=s c a l e L 1 / p i x l 1
p p i x l 2=s c a l e L 2 / p i x l 2
; f i n a l equivalent steps in p i x e l
; Ã¨ i l numero d i p i x e l e f f e t t i v o con i l q u a l e devo s p o s t a r e
N1step=Nstep L1 / p p i x l 1
N2step=Nstep L2 / p p i x l 2
i l frame l a y e r su o g n i p i x e l d e l l
k imm=1
n imm=200
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
Z d i m en s i o n =200
Z v e c t= f l t a r r ( Z d i m en s i o n )
; print , Z vect
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
N 1=12
; number o f subap i n t h e SH f i r s t an n u l u s
; number o f subap i n t h e SH s e c on d an n u l u s
N 2=18
; t o t a l number o f s u b a p e r t u r e i n t h e SH ar r ay
N Subap=N 1+N 2
Xm=60
; number o f columns i n each s u b a p e r t u r e
Ym=60
; number o f rows i n each s u b a p e r t u r e
Dim X=712
; number o f rows o f t h e SH ar r ay
Dim Y=632
Centr POS=i n t a r r ( 2 , N Subap ) ; s u b a p e r t u r e c e n t r e p o s i t i o n s i n t h e SH ar r ay
Centr PUP=i n t a r r ( 2 , N Subap ) ; s u b a p e r t u r e c e n t r e p o s i t i o n s i n t h e p u p i l ar r ay
Xb=64
; number o f columns i n s u b a p e r t u r e u s e d f o r c o r r e l a t i o n
Yb=64
; number o f rows i n s u b a p e r t u r e u s e d f o r c o r r e l a t i o n
x c o r r =16
; number o f columns i n c o r r e l a t i o n ar r ay
y c o r r =16
; number o f rows i n c o r r e l a t i o n ar r ay
MAX R=4
;####################################################
; ∗∗∗∗∗∗∗∗∗ ∗ ∗ ∗∗ ∗ ∗∗ ∗ ∗ ∗∗ ∗ ∗ ∗∗ ∗ ∗∗ ∗ ∗ ∗∗ ∗ ∗ ∗∗ ∗ ∗∗ ∗ ∗ ∗+
; 0. 475
; T e l e s c o p e a p e r t u r e w i d t h (m)
telescope D = 0.90
t e l e s c o p e O = 0 . 4 0 ; 0. 396
;0.0
; T e l e s c o p e c e n t r a l obscurament w i d t h (m)
lambda
= 5 0 0 0 . e −10
; 4305. e−10
; 4 5 0 7 . e−10
; Working w a v e l e n g t h (m)
pix
=
0.215
; P i x e l −s i z e en a r c s e c ( 0 . 0 8 2 9 3 7 a t 450 nm) ( 0 . 0 8 2 8 5 2 a t 430 nm)
104
s i z =500 ; 475
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
max=f i x ( n p i x / N1step )
N1step=f i x ( N1step )
N2step=f i x ( N2step )
print , max−1 , ’ time s t e p s ’
TOT Z VECT= f l t a r r ( Z d i m en s i o n , max )
SQR Z VECT= f l t a r r ( Z d i m en s i o n , max )
T o t v e c t= f l t a r r (2 ∗ N subap , max)
print , s y s t i m e ( 0 )
r e s t o r e , ’ / media / d i s k /G200 . s a v ’ , / v e r b o s e
G0=G( 0 : Z d i m en s i o n −1 ,∗)
G=G0
svdc ,G,W, U,V
; ı̈¿ 21 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21
y1=7+ f i x (ym / 2 . ) ; Cen x −2∗Xm−2∗ d i s t /2
y2=99+ f i x (ym/ 2 . ) ; Cen x −1∗Xm−1∗ d i s t /2
y3=192+ f i x (ym / 2 . ) ; Cen x
y4=284+ f i x (ym / 2 . ) ; Cen x+1∗Xm+1∗ d i s t /2
y5=376+ f i x (ym / 2 . ) ; Cen x+2∗Xm+2∗ d i s t /2
y6=468+ f i x (ym / 2 . ) ; Cen x+2∗Xm+2∗ d i s t /2
y7=559+ f i x (ym / 2 . ) ; Cen x+2∗Xm+2∗ d i s t /2
; ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21
x1=10+ f i x (xm/ 2 . ) ; Cen y −2∗Ym−2∗ d i s t
x2=62+ f i x (xm/ 2 . ) ; Cen y −1.5∗Ym−1∗ d i s t
x3=115+ f i x (xm / 2 . ) ; Cen y −1∗Ym−1∗ d i s t
x4=168+ f i x (xm / 2 . ) ; Cen y −0.5∗Ym−0.5∗ d i s t
x5=222+ f i x (xm / 2 . ) ; Cen y
x6=274+ f i x (xm / 2 . ) ; Cen y +0.5∗Ym+0.5∗ d i s t
x7=327+ f i x (xm / 2 . ) ; Cen y+1∗Ym+1∗ d i s t
x8=380+ f i x (xm / 2 . ) ; Cen y +1.5∗Ym+1∗ d i s t
x9=434+ f i x (xm / 2 . ) ; Cen y+2∗Ym+2∗ d i s t
x10=486+ f i x (xm / 2 . ) ; Cen y+2∗Ym+2∗ d i s t
x11=539+ f i x (xm / 2 . ) ; Cen y+2∗Ym+2∗ d i s t
x12=592+ f i x (xm / 2 . ) ; Cen y+2∗Ym+2∗ d i s t
x13=646+ f i x (xm / 2 . ) ; Cen y+2∗Ym+2∗ d i s t
; ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21
;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
;++++++++++++++++++++++++++++++++++++
; ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 FIRST ANNULUS
Centr POS ( ∗ , 1 8 ) = [ x7 , y6 ]
Centr POS ( ∗ , 1 9 ) = [ x9 , y6 ]
Centr POS ( ∗ , 2 0 ) = [ x10 , y5 ]
105
Centr POS ( ∗ , 2 1 ) = [ x11 , y4 ]
Centr POS ( ∗ , 2 2 ) = [ x10 , y3 ]
Centr POS ( ∗ , 2 3 ) = [ x9 , y2 ]
Centr POS ( ∗ , 2 4 ) = [ x7 , y2 ]
Centr POS ( ∗ , 2 5 ) = [ x5 , y2 ]
Centr POS ( ∗ , 2 6 ) = [ x4 , y3 ]
Centr POS ( ∗ , 2 7 ) = [ x3 , y4 ]
Centr POS ( ∗ , 2 8 ) = [ x4 , y5 ]
Centr POS ( ∗ , 2 9 ) = [ x5 , y6 ]
; ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 SECOND ANNULUS
Centr POS ( ∗ , 0 ) = [ x8 , y7 ]
Centr POS ( ∗ , 1 ) = [ x10 , y7 ]
Centr POS ( ∗ , 2 ) = [ x11 , y6 ]
Centr POS ( ∗ , 3 ) = [ x12 , y5 ]
Centr POS ( ∗ , 4 ) = [ x13 , y4 ]
Centr POS ( ∗ , 5 ) = [ x12 , y3 ]
Centr POS ( ∗ , 6 ) = [ x11 , y2 ]
Centr POS ( ∗ , 7 ) = [ x10 , y1 ]
Centr POS ( ∗ , 8 ) = [ x8 , y1 ]
Centr POS ( ∗ , 9 ) = [ x6 , y1 ]
Centr POS ( ∗ , 1 0 ) = [ x4 , y1 ]
Centr POS ( ∗ , 1 1 ) = [ x3 , y2 ]
Centr POS ( ∗ , 1 2 ) = [ x2 , y3 ]
Centr POS ( ∗ , 1 3 ) = [ x1 , y4 ]
Centr POS ( ∗ , 1 4 ) = [ x2 , y5 ]
Centr POS ( ∗ , 1 5 ) = [ x3 , y6 ]
Centr POS ( ∗ , 1 6 ) = [ x4 , y7 ]
Centr POS ( ∗ , 1 7 ) = [ x6 , y7 ]
; ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12
;++++++++++++++++++++++++++++++++++++++++++++++++++++++++
f o r mm=0 ,200 do begin
p h a s e r e f =(s cr een G ( 0 : 1 0 4 3 , ( N1step ∗mm) : ( N1step ∗mm)+1043)+ s c r e e n L ( (mm∗ N2step ) :mm∗ N2step + 1 0 4 3 , 0 :
REFERENCE= f l t a r r ( Dim X , Dim Y )
imm=p h a s e r e f
For i s u b a p =0 , N subap −1 do begin
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
xpos=Centr POS ( 0 , i s u b a p )− f i x (xm / 2 . )
ypos=Centr POS ( 1 , i s u b a p )− f i x (ym / 2 . )
sub imm=imm( xpos : xpos+xm−1 , ypos : ypos+ym−1)
REFERENCE( xpos : xpos+xm−1 , ypos : ypos+ym−1) =sub imm
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
endfor
;;####################################################
; window , 0 , x s i z e=Dim X , y s i z e=Dim Y , t i t l e =’REFERENCE’
; t v s c l , imm;−REFERENCE
;;####################################################
106
vx=i n t a r r (xm)
vy=i n t a r r (ym)
For x=0 ,xm−1 do vx ( x)=x−xm/2
For y=0 ,ym−1 do vy ( y)=y−ym/2
vx= s h i f t ( vx , xm/ 2 )
vy= s h i f t ( vy , ym/ 2 )
Mx=vx#r e p l i c a t e ( 1 ,ym)
My=r e p l i c a t e ( 1 ,xm)#vy
For x=0 ,xm−1 do vx ( x)=x
For y=0 ,ym−1 do vy ( y)=y
Px=vx#r e p l i c a t e ( 1 ,ym)
Py=r e p l i c a t e ( 1 ,xm)#vy
v s h i f t= f l t a r r ( 2 , N Subap )
; ˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆ
; Fr e c u e n c i a de c o r t e en radˆ−1
n u c u t o f f=t e l e s c o p e D / lambda
pixrad = pix ∗! pi /180./3600.
; P i x e l −s i z e en r a d i a n .
d e l t a n u =1./ s i z / p i x r a d
; I n t e r v a l o de muestreo en radˆ−1
r p u p i l=n u c u t o f f / d e l t a n u / 2 .
; Radio de l a p u p i l a en p i x e l s
p i x c u t =2∗ r p u p i l
s i z =2.ˆ f i x ( a l o g 1 0 ( p i x c u t )/ a l o g 1 0 ( 2 . ) +1) +20 ; l a prima m a t r i c e con dim 2ˆn che c o n t i e n e l a p u p i l l a
s h i f t x= f l t a r r ( N Subap )
s h i f t y= f l t a r r ( N Subap )
ARRAY= f l t a r r ( Dim x , Dim y )
PSF ARRAY= f l t a r r ( Dim x , Dim y )
PSF ref ARRAY= f l t a r r ( Dim x , Dim y )
D v ect= f l t a r r (2 ∗ N subap )
; rd=radius aper ANN ( s i z /2 , p i x c u t , t e l e s c o p e O / t e l e s c o p e D , t s u p p o r t )
; ˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆ;ˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆ
;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
Cen=f i x ( s i z / 2 . )
e x a g s i d e= 84 ; 7 9 ; 1 0 9 ; 2 0 0
e x a g a p t=e x a g s i d e ∗ c o s ( ! PI / 6 . )
; ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21
y1=Cen −4.5∗ e x a g s i d e
y2=Cen−3∗ e x a g s i d e
y3=Cen −1.5∗ e x a g s i d e
y4=Cen
y5=Cen +1.5∗ e x a g s i d e
y6=Cen+3∗ e x a g s i d e
y7=Cen +4.5∗ e x a g s i d e
; ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21
x1=Cen−6∗ e x a g a p t
x2=Cen−5∗ e x a g a p t
x3=Cen−4∗ e x a g a p t
x4=Cen−3∗ e x a g a p t
x5=Cen−2∗ e x a g a p t
x6=Cen−1∗ e x a g a p t
x7=Cen
x8=Cen+1∗ e x a g a p t
x9=Cen+2∗ e x a g a p t
x10=Cen+3∗ e x a g a p t
x11=Cen+4∗ e x a g a p t
x12=Cen+5∗ e x a g a p t
x13=Cen+6∗ e x a g a p t
; ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21
;@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
;++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 FIRST ANNULUS
Centr PUP ( ∗ , 1 8 ) = [ x7 , y6 ]
107
Centr PUP ( ∗ , 1 9 ) = [ x9 , y6 ]
Centr PUP ( ∗ , 2 0 ) = [ x10 , y5 ]
Centr PUP ( ∗ , 2 1 ) = [ x11 , y4 ]
Centr PUP ( ∗ , 2 2 ) = [ x10 , y3 ]
Centr PUP ( ∗ , 2 3 ) = [ x9 , y2 ]
Centr PUP ( ∗ , 2 4 ) = [ x7 , y2 ]
Centr PUP ( ∗ , 2 5 ) = [ x5 , y2 ]
Centr PUP ( ∗ , 2 6 ) = [ x4 , y3 ]
Centr PUP ( ∗ , 2 7 ) = [ x3 , y4 ]
Centr PUP ( ∗ , 2 8 ) = [ x4 , y5 ]
Centr PUP ( ∗ , 2 9 ) = [ x5 , y6 ]
; ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 SECOND ANNULUS
Centr PUP ( ∗ , 0 ) = [ x8 , y7 ]
Centr PUP ( ∗ , 1 ) = [ x10 , y7 ]
Centr PUP ( ∗ , 2 ) = [ x11 , y6 ]
Centr PUP ( ∗ , 3 ) = [ x12 , y5 ]
Centr PUP ( ∗ , 4 ) = [ x13 , y4 ]
Centr PUP ( ∗ , 5 ) = [ x12 , y3 ]
Centr PUP ( ∗ , 6 ) = [ x11 , y2 ]
Centr PUP ( ∗ , 7 ) = [ x10 , y1 ]
Centr PUP ( ∗ , 8 ) = [ x8 , y1 ]
Centr PUP ( ∗ , 9 ) = [ x6 , y1 ]
Centr PUP ( ∗ , 1 0 ) = [ x4 , y1 ]
Centr PUP ( ∗ , 1 1 ) = [ x3 , y2 ]
Centr PUP ( ∗ , 1 2 ) = [ x2 , y3 ]
Centr PUP ( ∗ , 1 3 ) = [ x1 , y4 ]
Centr PUP ( ∗ , 1 4 ) = [ x2 , y5 ]
Centr PUP ( ∗ , 1 5 ) = [ x3 , y6 ]
Centr PUP ( ∗ , 1 6 ) = [ x4 , y7 ]
Centr PUP ( ∗ , 1 7 ) = [ x6 , y7 ]
; ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12 ı̈¿ 21 ı̈¿ 12 ı̈¿ 12
;++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Exag mask=i n t a r r ( e x a g s i d e ∗2+1 , e x a g s i d e ∗2+1)
window , 3 , x s i z e=e x a g s i d e ∗2+1 , y s i z e=e x a g s i d e ∗2+1 , t i t l e = ’ Exag mask ’
t v s c l , Exag mask
xpos=e x a g s i d e
ypos=e x a g s i d e
y1=Ypos−e x a g s i d e
y2=Ypos−e x a g s i d e /2
y3=Ypos+e x a g s i d e /2
y4=Ypos+e x a g s i d e
x1=Xpos−e x a g a p t
x2=Xpos
x3=Xpos+e x a g a p t
p l o t s , [ X2 , X3 ] , [ Y1 , Y2 ] , / d e v i c e , c o l o r = [ 2 5 5 , 2 5 5 , 2 5 5 ] , Thick=2
p l o t s , [ X3 , X3 ] , [ Y2 , Y3 ] , / d e v i c e , c o l o r = [ 2 5 5 , 2 5 5 , 2 5 5 ] , Thick=2
p l o t s , [ X3 , X2 ] , [ Y3 , Y4 ] , / d e v i c e , c o l o r = [ 2 5 5 , 2 5 5 , 2 5 5 ] , Thick=2
p l o t s , [ X2 , X1 ] , [ Y4 , Y3 ] , / d e v i c e , c o l o r = [ 2 5 5 , 2 5 5 , 2 5 5 ] , Thick=2
p l o t s , [ X1 , X1 ] , [ Y3 , Y2 ] , / d e v i c e , c o l o r = [ 2 5 5 , 2 5 5 , 2 5 5 ] , Thick=2
p l o t s , [ X1 , X2 ] , [ Y2 , Y1 ] , / d e v i c e , c o l o r = [ 2 5 5 , 2 5 5 , 2 5 5 ] , Thick=2
b=LABEL REGION(255 − t v r d ( ) , /EIGHT)
Exag mask=b eq 3
w d el ete , 3
; ˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆ;ˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆˆ
;FOR Z J =2, Z dimension −1 DO BEGIN
Z v e c t=Z v e c t ∗0
REFERENCE= f l t a r r ( Dim X , Dim Y )
imm=p h a s e r e f
108
For i s u b a p =0 , N subap −1 do begin
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
xpos=Centr POS ( 0 , i s u b a p )− f i x (xm / 2 . )
ypos=Centr POS ( 1 , i s u b a p )− f i x (ym / 2 . )
sub imm=imm( xpos : xpos+xm−1 , ypos : ypos+ym−1)
REFERENCE( xpos : xpos+xm−1 , ypos : ypos+ym−1) =sub imm
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
endfor
;;####################################################
; window , 0 , x s i z e=Dim X , y s i z e=Dim Y , t i t l e =’REFERENCE’
; t v s c l ,REFERENCE
;
;;####################################################
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
d e v i c e , decomposed=0
l o a d c t , 3 , / SILENT
window , 1 , x s i z e=s i z , y s i z e=s i z , t i t l e = ’PHASE ’ , xpos =0 , Ypos=0
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
phase=p h a s e r e f
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
wset , 1
t v s c l , phase
phase=complex ( c o s ( phase ) , s i n ( phase ) )
arrow , 5 0 , 0 , 5 0 , 2 ∗mm
arrow , 0 , 5 0 ,mm, 5 0
wset , 1
;DISEGNO L ’ARRAY DI SOTTOPUPILLE
For I s u b a p =0 , N subap− 1 do begin
xpos=Centr PUP ( 0 , I s u b a p )
ypos=Centr PUP ( 1 , I s u b a p )
y1=Ypos−e x a g s i d e
y2=Ypos−e x a g s i d e /2
y3=Ypos+e x a g s i d e /2
y4=Ypos+e x a g s i d e
x1=Xpos−e x a g a p t
x2=Xpos
x3=Xpos+e x a g a p t
p l o t s , [ X2 , X3 ] , [ Y1 , Y2 ] , / d e v i c e , c o l o r =255 ; [ 0 , 0 , 0 ]
p l o t s , [ X3 , X3 ] , [ Y2 , Y3 ] , / d e v i c e , c o l o r =255 ; [ 0 , 0 , 0 ]
p l o t s , [ X3 , X2 ] , [ Y3 , Y4 ] , / d e v i c e , c o l o r =255 ; [ 0 , 0 , 0 ]
p l o t s , [ X2 , X1 ] , [ Y4 , Y3 ] , / d e v i c e , c o l o r =255 ; [ 0 , 0 , 0 ]
p l o t s , [ X1 , X1 ] , [ Y3 , Y2 ] , / d e v i c e , c o l o r =255 ; [ 0 , 0 , 0 ]
p l o t s , [ X1 , X2 ] , [ Y2 , Y1 ] , / d e v i c e , c o l o r =255 ; [ 0 , 0 , 0 ]
p l o t s , [ Xpos , Xpos ] , [ Ypos+5 ,Ypos − 5 ] , / d e v i c e , c o l o r =255 ; [ 0 , 0 , 0 ]
p l o t s , [ Xpos −5 ,Xpos + 5 ] , [ Ypos , Ypos ] , / d e v i c e , c o l o r =255 ; [ 0 , 0 , 0 ]
endfor
i d =0
; N subap− 1
For I s u b a p =0 , N subap −1 do begin
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
109
xpos=Centr POS ( 0 , i s u b a p )− f i x (xm / 2 . )
ypos=Centr POS ( 1 , i s u b a p )− f i x (ym / 2 . )
imm=REFERENCE( xpos : xpos+xm−1 , ypos : ypos+ym−1)
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
xpos=Centr PUP ( 0 , I s u b a p )
ypos=Centr PUP ( 1 , I s u b a p )
subpup= f l t a r r ( e x a g s i d e ∗ 2 , e x a g s i d e ∗ 2 )
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
MASK=Exag mask ∗SUPPORT( xpos−e x a g s i d e : xpos+e x a g s i d e
;
;
subpup=MASK∗PHASE( xpos−e x a g s i d e
window , 3
t v s c l , subpup
: xpos+e x a g s i d e
, ypos−e x a g s i d e : ypos+e x a g s i d e )
, ypos−e x a g s i d e : ypos+e x a g s i d e )
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
F IMM= f f t (imm, −1)
r= f f t (MASK, −1)
c r o s s= f f t ( c o n j ( r ) ∗ r , 1 )
c r o s s=c r o s s / c r o s s ( 0 , 0 )
P S F r ef=f l o a t ( f f t ( c r o s s , 1 ) )
P S F r ef=s h i f t ( PSF ref , e x a g s i d e , e x a g s i d e )
; c o n v o l u t i o n i n f o u r i e r s pac e
r e f= f f t (F IMM∗ c r o s s , 1 )
f= f f t ( subpup , −1)
c r o s s= f f t ( c o n j ( f ) ∗ f , 1 )
c r o s s=c r o s s / c r o s s ( 0 , 0 )
PSF=f l o a t ( f f t ( c r o s s , 1 ) )
PSF=s h i f t (PSF , e x a g s i d e , e x a g s i d e )
subap= f f t (F IMM∗ c r o s s , 1 )
; r e s u l t Ã¨ l a p s f t e o r i c a f i t t a t a con una g a u s s i a n a
R e s u l t = GAUSS2DFIT (PSF , A1 )
R e s u l t = GAUSS2DFIT ( PSF ref , A2 )
s h i f t x ( I s u b a p )=(A1(4) −A2 ( 4 ) )
s h i f t y ( I s u b a p )=(A1(5) −A2 ( 5 ) )
; print , s h i f t x ( I subap ) , s h i f t y ( I subap )
D v ect ( i d )= s h i f t x ( i s u b a p )
D v ect ( i d +1)= s h i f t y ( i s u b a p )
i d=i d +2
i f abs ( s h i f t x ( i s u b a p ) ) gt 0 . 0 or abs ( s h i f t y ( i s u b a p ) ) gt 0 . 0 then begin
t a i l =[ xpos , ypos ]
h ea d =[ xpos +N1step ∗ s h i f t x ( I s u b a p ) , ypos +N2step ∗ s h i f t y ( I s u b a p ) ]
wset , 1
ARROW, t a i l ( 0 ) , t a i l ( 1 ) , h ea d ( 0 ) , h ea d ( 1 ) , c o l o r =255 , / s o l i d
, h t h i c k = 2 . 0 , t h i c k =2.
endif
endfor
T o t v e c t ( ∗ ,mm)= D v ect
Z v e c t= SVSOL(U, W, V, D v ect )
Z vect (1)=0.
110
TOT Z VECT( ∗ ,mm)= Z v e c t / ( 2 ∗ ! p i )
z v e c t=z v e c t / ( 2 ∗ ! p i )
stop
w d el ete , 1
d v e c t=d v e c t ∗0
z v e c t=z v e c t ∗0
mask=mask ∗0
exag mask=exag mask ∗0
s i z =500
endfor
stop
t o t z v e c t t s=t o t z v e c t
save , f i l e n a m e= ’ / media / d i s k / t o t z v e c t t s 2 . s a v ’ , t o t z v e c t t s
save , f i l e n a m e= ’ / media / d i s k / s h i f t s . s a v ’ , t o t v e c t
end
pro PHASE SCREEN
; This program g e n e r a t e s l a r g e phas e s c r e e n w i t h Kolmogorov or Von−Karman t u r b u l e n c e power spectrum
;
;
;
;
;
;
;
;
− L0 i s t h e o u t e r s c a l e where e n e r g y i s i n j e c t e d ( i t i s i n f i n i t y f o r Kolmogorov model )
− r0 Fr i e d parameter
− Nx , Ny ar e t h e di m e n s i on s o f phas e s c r e e n s u p p o r t
− Gx , Gy ar e t h e di m e n s i on s o f phas e s c r e e n i n m e t e r s
− D=G/N i s t h e s am pl i n g r a t i o ( eg . 1 p i x e l=D m e t e r s )
− spectrum =1,0
; SPECTRUM=1 VON−KARMAN
SPECTRUM=0 KOLMOGOROV
− ” s c r e e n ” i s t h e f i n a l ar r ay c o n t a i n i n g s i m u l a t e d phas e
Marco S t a n g a l i n i 2008
r e s t o r e , ’ / media / u s b d i s k / s u p p o r t . s a v ’
DEVICE, DECOMPOSED = 0
loadct ,13
Nx=80
Ny=80
r 0 =0.15
; Outer s c a l e i n m e t e r s
Gx=1.8
Gy=1.8
; Fr i e d parameter i n m e t e r s
; Phase s c r e e n di m e n s i on s i n m e t e r s
; Phase s c r e e n di m e n s i on s i n m e t e r s
; ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ TURBULENCE SPECTRUM TYPE ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
111
spectrum =1
L0 =.01
; SPECTRUM=1 VON−KARMAN
SPECTRUM=0 KOLMOGOROV
i f spectrum eq 1 then begin
s c a l e = 1 . / ( ( L0 ) ˆ ( 2 . ) )
print , ’ Von−Karman Spectrum ’
endif e l s e begin
s c a l e =0.
print , ’ Kolmogorov Spectrum ’
endelse
print , ’ s p a t i a l s a m p l i n g : ’ , Gx/Nx
ydim=Ny
xdim=Nx
s e e d=SYSTIME ( 1 )
r e a l =(RANDOMN( s eed , ydim , xdim , / normal ) )
; im=(RANDOMN( seed , ydim , xdim , / normal ) )
H= f f t ( r e a l , −1)
P h a s e r e= f l t a r r (Ny , Nx)
Phase im= f l t a r r (Ny , Nx)
; Phase=complex ( P has e r e , Phase im )
phase=co m p l ex a r r ( ny , nx )
F= f l t a r r ( ny , nx )
T=i n t a r r ( ny , nx )
phase=phase ∗0
p h a s e r e=p h a s e r e ∗ 0 .
p h a s e i m=p h a s e i m ∗ 0 .
; ANTIALIASING FILTER
For t1 =0 ,Ny−1 do begin
For z1 =0 ,Nx−1 do begin
i f ( ( ( z1−Nx / 2 . ) / ( Nx / 2 . ) ) ˆ ( 2 . ) + ( ( t1−Ny / 2 . ) / ( Ny / 2 . ) ) ˆ ( 2 . ) ) l e 1 then begin
T( t1 , z1 )=1
endif e l s e begin
T( t1 , z1 )=0
endelse
Endfor
Endfor
; Power spectrum Von Karman
f o r m1=0 ,Nx−1 do begin
f o r n1 =0 ,Ny−1 do begin
i f ( (m1−Nx/ 2 ) eq 0 ) AND ( ( n1−Ny/ 2 ) eq 0 ) then begin
F( n1 , m1)=0
endif e l s e begin
F( n1 , m1)= 0 . 1 5 1 3 ∗ ( Gx∗Gy ) ˆ ( 1 . / 2 . ) ∗ r 0 ˆ ( − 5 . / 6 . ) ∗ ( ( (m1−Nx/ 2 )/ Gx ) ˆ ( 2 . ) + ( ( n1−Ny/ 2 )/Gy ) ˆ ( 2 . )
endelse
endfor
endfor
;
f=g a u s s 2 d f i t (F)
; r e a l p a r t Noise&VonKarman c o n v o l u t i o n
p h a s e r e=H∗F∗T
r e=t o t a l ( p h a s e r e )
FOR m=0 ,Nx−1
DO BEGIN
112
print , ’ countdown ’ , Nx−m
FOR n=0 ,Ny−1 DO BEGIN
f o r m1=0 ,Nx−1 do begin
f o r n1 =0 ,Ny−1 do begin
p h a s e i m ( n ,m) = ( 2 . ∗ ! p i ∗ ( ( ( m1−Nx / 2 ) ∗ (m−Nx / 2 ) ) / Nx+(( n1−Ny / 2 ) ∗ ( n−Ny / 2 ) ) / Ny ) ) ∗ i m a g i n a r y (H( n ,m) )
endfor
endfor
phase ( n ,m)=complex ( r e ∗ c o s ( t o t a l ( Phase im ) ) , ( r e ∗ s i n ( t o t a l ( p h a s e i m ) ) ) )
ENDFOR
ENDFOR
s c r e e n =( f f t ( phase , 1 ) )
window , 0 , t i t l e = ’SIMULATED PHASE SCREEN ’ ; , x s i z e =700, y s i z e =680
s h a d e s u r f , s c r e e n , x c h a r s i z e =2 , y c h a r s i z e =2 , z c h a r s i z e=2
; ;∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
;;
; Ns=3
; subharmonic k e r n e l s i z e
; Np=19
; subharmonic number
; r e a l 2 =(RANDOMN( seed , 2 ∗ Ns , 2 ∗ Ns , Np, / normal ) )
; a=0.
; const
; wp=1.
; weight f ac t or
;
; Fsub=f l t a r r ( 2∗ Ns , 2 ∗ Ns , Np)
; Hsub=f l t a r r ( 2∗ Ns , 2 ∗ Ns , Np)
;
; f o r p=0,Np−1 do b e g i n
; Hsub ( ∗ , ∗ , p)= f f t ( r e a l 2 ( ∗ , ∗ , p ) , −1)
; endfor
;
; p=0
;
;
f o r p=0,Np−1 do b e g i n
;
f o r m1=0,Ns−1 do b e g i n
;
f o r n1=0,Ns−1 do b e g i n
;
i f ( (m1−Ns /2) eq 0) AND ( ( n1−Ns /2) eq 0) t he n b e g i n
;
Fsub ( n1 , m1, p)=0
;
endif e lse begin
;
Fsub ( n1 , m1, p )=0.1513∗( Gx∗Gy ) ˆ ( 1 . / 2 . ) ∗ r0 ˆ( − 5. /6. ) ∗ 3. ˆ( − p ) ∗ ( ( 3. ˆ( − p ) ∗ ( ( m1−Ns/2)+a )/Gx) ˆ( 2. ) + ( 3.
;
endelse
;
endfor
;
endfor
;
Fsub ( ∗ , ∗ , p)= g a u s s 2 d f i t ( Fsub ( ∗ , ∗ , p ) )
;
endfor
;
;
;
;
; p=0
; p h a s e r e s u b= f l t a r r ( 2∗ Ns , 2 ∗ Ns )
; p h a s e i m s u b= f l t a r r (Ny , Nx , Np)
; p h a s e s u b= f l t a r r (Ny , Nx)
;
;
f o r p=0,Np−1 do b e g i n
;
p h a s e r e s u b (∗ ,∗)= p h a s e r e s u b (∗ ,∗)+ Fsub ( ∗ , ∗ , p ) ∗ Hsub ( ∗ , ∗ , p )
;
endfor
;
; b=f i x (Ny/Ns )
; r e s u b=t o t a l ( p h a s e r e s u b )
113
;
; For p=0,Np−1 do b e g i n
;FOR r =0,Nx−1, b DO BEGIN
; FOR s =0,Ny−1, b DO BEGIN
;
f o r m1=0,Ns−1 do b e g i n
;
f o r n1=0,Ns−1 do b e g i n
;
p h a s e i m s u b ( r , s , p ) = ( 2 . ∗ ! p i ∗3.ˆ( − p ) ∗ ( ( ( m1+(Ns/2)+a )∗ r /Nx)+(( n1+(Ns/2)+a )∗ s /Ny ) ) ) ; ∗ i m
;
endfor
;
endfor
;
p h a s e s u b ( r , s )= p h a s e s u b ( r , s )+complex ( r e s u b ∗ c os ( t o t a l ( P has e i m s u b ( ∗ , ∗ , p ) ) ) , ( r e s u b ∗ s
; ENDFOR
;ENDFOR
; endfor
;
; p=0
; phas e s u b har m on i c=p h a s e s u b ( ∗ , ∗ ) ∗ 0
;
;
;
; i n t= f f t (T∗ p h a s e s u b , 1 )
; ; i n t=i n t e r p o l a t e ( phas e s u b har m on i c , where ( phas e s u b har m on i c l t 1 0 . ∗ max ( phas e s u b har m on i c ) ) )
; window , 2
; shade surf , int ∗ screen
window , 0 , x s i z e =600 , y s i z e =600
expand , s c r e e n ( 5 : 7 5 , 5 : 7 5 ) , 2 0 8 8 , 2 0 8 8 , s c r e e n
tvs cl , congrid ( screen ,600 ,600)
; save , f i l e n a m e =’/media / u s b d i s k / s c r e e n L . sav ’ , s c r e e n
stop
end
pro t i p t i l t h i g h O r d e r
; This program compute h i g h o r d e r Z e r n i k e modes c o n t r i b u t i o n t o Tip T i l t
DEVICE, DECOMPOSED = 0
loadct , 3
r e s t o r e , ’ / media/ u s b d i s k / z e r n i k e 2 0 0 G f l o a t . s a v ’ , / v e r b o s e
r e s t o r e , ’ / media/ u s b d i s k / t o t z v e c t 2 0 0 . s a v ’ , / v e r b o s e
p h a s e t i p t i l t= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
phase HO= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
p h a s e r e c= f l t a r r ( 1 0 4 4 , 1 0 4 4 )
r e c z e= f l t a r r ( 2 0 0 )
immagine=2
;USAGE t o t z v e c t ( r , immagine )
f o r i =4 ,199 do begin
p h a s e t i p t i l t=p h a s e t i p t i l t+t o t z v e c t ( i , immagine ) ∗ z e r ( ∗ , ∗ , i )
print , i
endfor
114
phase HO= p h a s e t i p t i l t
; ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ PROIEZIONE WF ON ZE ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
f o r r =2 ,3 do begin
r e c z e ( r )= t o t a l ( ( ( ZER( ∗ , ∗ , r ) ) ∗ phase HO )/ t o t a l ( z e r ( ∗ , ∗ , 0 ) ) )
phase HO=phase HO−r e c z e ( r )∗ZER( ∗ , ∗ , r )
; p r i n t , t o t z v e c t ( r , immagine ) , r e c z e ( r )
endfor
; ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
f o r t =2 ,3 do begin
p h a s e r e c=p h a s e r e c+r e c z e ( t )∗ z e r ( ∗ , ∗ , t )
endfor
;
; SET PLOT, ’PS ’
;
; ; S e t t h e P o s t S c r i p t d e v i c e t o ∗8∗ b i t s pe r c o l o r , n ot 2 4 .
;
; DEVICE, FILE=’/home/ s o l a r −x2 / Desktop / t i p t i l t r e s i d u a l . ps ’ , /COLOR, BITS=8
;WINDOW, 0 , TITLE=’ R e f e r c e n c e ’
; shade surf , p h a s e t i p t i l t
; window , 1 , t i t l e =’ High o r d e r modes t i p − t i l t c o n t r i b u t i o n ’
s u r f a c e , p h a s e r e c , z t i t l e = ’ Phase l a g [ waves ] ’
; DEVICE, /CLOSE
; ; Return p l o t t i n g t o X windows :
; SET PLOT, ’X’
stop
end
pro t s a n a l i s y s
; Z e r n i k e f o r e c a s t i n g w i t h l i n e a r p r e d i c t i o n method b as e d upon m at r i x c a l c u l a t i o n
DEVICE, DECOMPOSED = 0
r e s t o r e , ’ / media / d i s k / t o t z v e c t 2 0 0 . s a v ’ , / v e r b o s e
; t o t z v e c t t s ( z e r n i k e , immagine )
l e a r n i n g =20
; t o t z v e c t=t o t z v e c t t s
l a s t= f l t a r r ( l e a r n i n g )
f o r ze =4 ,190 do begin
F= f l t a r r ( l e a r n i n g , l e a r n i n g )
H= f l t a r r ( l e a r n i n g , l e a r n i n g )
; f u t u r e m at r i x
; h i s t o r i c a l m at r i x
f o r t =0 , l e a r n i n g −1 do begin
F( t ,∗)= t o t z v e c t ( ze , t +1)
115
H( t ,∗)= t o t z v e c t ( ze , t )
endfor
; P=F#t r a n s p o s e ( c o n j (H) )
P=(F#t r a n s p o s e (H))# (( i n v e r t (H#t r a n s p o s e (H ) ) ) )
f o r l =0 , l e a r n i n g −1 do begin
l a s t ( l )=H( 0 , l )
endfor
s=r e a l p a r t (P#l a s t )
print , r e a l p a r t ( s (0)+ s ( 1 ) )
print , t o t z v e c t ( ze , l e a r n i n g +1)
endfor
;
; max ze=90
; x= f l t a r r ( 200)
; lagmax=90
; sh=f l t a r r ( max ze )
;H= f l t a r r ( max ze , max ze )
; F= f l t a r r ( max ze , max ze )
; l e a r n i n g =50
;
;
;
;
f o r t =0,89 do b e g i n
;
x ( t )= t o t z v e c t t s ( ze , t )
;
endfor
;
; t i m e s e r i e s= f l t a r r ( 220)
;ACF=IMSL AUTOCORRELATION( x , lagmax )
;
; window , 1 , t i t l e =’SAMPLE AUTOCOVARIANCE’
; p l o t ,ACF, y t i t l e =’ACF’ , x t i t l e =’Lag [ ms ] ’ , psym =5; , x r an g e = [ 0 : 5 0 ]
;
; arma = IMSL ARMA( x , 1 0 , 1 , n p r e d i c t =50, f o r e c a s t=f o r e )
;
;
;
f o r t =0, max ze −1 do b e g i n
;
sh ( t )= t o t z v e c t t s ( ze , t )
;
endfor
;
;
; P=(F#t r a n s p o s e (H))#(( i n v e r t (H#t r a n s p o s e (H) ) ) )
;
; ! P .MULTI = [ 0 , 1 , 2 ]
; window , 6 , t i t l e =’ For e c as t ’ , x s i z e =1000
; p l o t , t o t z v e c t t s ( 0 : max ze −1 ,97)∗10
; loadct ,3
116
; o p l o t , ( P##( t o t z v e c t t s ( 0 : max ze − 1 , 96) ) ) , c o l o r =200; , psym=2
; loadct ,0
; p l o t , ( ( t o t z v e c t t s ( 0 : max ze −1,97) −(P##( t o t z v e c t t s ( 0 : max ze − 1 , 9 6 ) ) ) ) ˆ ( 2 . ) ) ˆ ( 0 . 5 ) , y r an g e = [ 0 , 1 ]
;
stop
end
117
118
Ringraziamenti
Desidero ringraziare il Prof. Francesco Berrilli e il Dott. Dario Del Moro per il loro supporto e la fiducia che
mi hanno sempre dimostrato già prima dell’inizio di questo lavoro di tesi. Una fiducia che mi ha stimolato
a fare sempre meglio. Voglio anche ringraziare gli amici del gruppo di Fisica Solare (Bartolomeo, Roberto,
Silvia e Valentina) per le chiacchierate e le risate che mi hanno fatto fare in questi mesi, Arnaldo che con i
suoi cioccolatini ha contribuito ad addolcire il lavoro, la mia famiglia e tutti i miei amici che hanno saputo
accompagnarmi attraverso i giorni più difficili.
Un sentito ringraziamento anche al Dott. Giuseppe Consolini per le stimolanti chiacchierate da cui ho
potuto trarre numerosi insegnamenti e spunti.
119
120
Pubblicazioni
121
Mem. S.A.It. Vol. 75, 282
c SAIt
2008
Memorie
della
Wavefront sensing for Themis: a case study for next
generation solar telescopes like EST
M. Stangalini,1 D. Del Moro,1 F. Berrilli1
Università di Roma ’Tor Vergata’ – Dipartimento di Fisica, Via della Ricerca
Scientifica 1, 00133 Rome, Italy
e-mail: [email protected]
Abstract. Wavefront sensing is one of the most challenging problem in solar
Adaptive Optics (AO) due to the presence of extended and low contrast AOtargets on the solar scene. We show preliminary results of wavefront modal phase
estimation on THEMIS Solar Telescope data using two approches: Zernike (ZE)
and Karhunen-Loeve (K-L) polynomial expansions. We have analyzed ShackHartmann (SH) dataset, and the wavefront phases were recontructed. A study
of pros and cons of the two different modal expansions has been carried out.
These results could be particularly helpful in developing and studying design and
performances for next generation solar telescopes based upon annular pupils.
Key words. Instrumentation: adaptive optics – Techniques: high angular resolution
1. Introduction
Solar AO, in contrast to non-solar case, has
to work on extended and poorly contrasted
targets.
In order to achieve a good wavefront description correlation technique has to be
carried out on SH subpupils images in order to estimate local wavefront slope from
their respective image shifts.
We present preliminary results of wavefront
sensing at Themis Solar Telescope using a
SH sensor.
Fig. 1. SH setup (30 subpupils) with measured
shifts superimposed (Black arrows).
Send offprint requests to: M. Stangalini
M. Stangalini et al.: Wavefront sensing for Themis
283
Fig. 2. (Left) Wavefront reconstruction with up to 200 ZE modes. (Centre) The same phase
reconstructed with up to 10 ZE polynomials. (Right) The same wavefront reconstructed with
only 10 K-L modes. It is worth noting that the first K-L mode contains already most of the
high spatial frequency features. Tip and tilt components have been forced to be ZE modes for
convenience in K-L expansion.
2. Dataset and wavefront slope
measurements
In this work we used SH sensor data acquired at Themis Solar Telescope on 3
August 2006 and slope measurements have
been estimated through Sum of Absolute
Differences correlation method on 30 SH
subpupils (Fig.1)
3. Modal wavefront reconstruction
To describe accurately wavefront aberrations we used a modal phase reconstruction based, at first, upon ZE polynomials (Noll, R.J. 1976) and then on KL functions (G.M. Dai 1996), both defined over annular pupil. The study of this
particular optical scheme is helpful in developing next generation European Solar
Telescope, giving us the possibility to test
wavefront reconstruction performance on
annular pupils with large central obscuration, tipical of Ritchey Chretien optical
configuration.
K-L decomposition makes use of Principal
Component Analisys to choose a suitable
subset of eigenvectors constructed as a
combination of ZE modes.
K-L, in contrast to ZE polynomials, are
statistically indipendent and so they represent a better basis for phase reconstruc-
tion. Fig.2 shows that using the same number of components, K-L functions (right)
results in a quite better wavefront description with respect to ZE polynomials (centre). This method ensure very high performance, speeding up realtime calculation
while maintaining the same accuracy.
4. Recostruction fitting error
Given a reference phase, built with up to
200 ZE mode, we tried to reconstruct this
phase using ZE and K-L approaches with
different number of modes. This analisys
has underlined how K-L modes are more
efficient in describing wavefront even at low
order modes with respect to ZE ones. At
low order K-L fitting error is always lower
than ZE reconstruction error.
References
Noll, R. J. (1976), Zernike polynomials
and atmospheric turbulence, J. Opt. Soc.
Am. 66, 207-211
G. M. Dai, ”Modal wave-front reconstruction with Zernike polynomials and
Karhunen-Loeve functions,” J. Opt. Soc.
Am. A 13, 1218- (1996)
12th European Solar Physics Meeting
Freiburg, Germany, 8 - 12 September 2008
Poster session 2.2
Improvements in Solar Adaptive Optics Correction using Short-time
Turbulence Forecasting
Stangalini, M.; Del Moro, D.; Berrilli, F.
Università di Roma Tor Vergata, ITALY
Abstract
Time delay error and integration time effects play a crucial role in Solar Adaptive Optics systems.
By using a multi-layer turbulence simulation we studied temporal variations of atmospheric turbulence.
We simulated Shack-Hartmann wavefront sensing with Zernike and Kharunen-Loeve modal
expansions too.
These two sets of polynomial basis have been defined over an annular pupil with a large central
obscurament, which is a common feature in Ritchey-Chretien optical configuration.
We found that Karhunen-Loeve functions, defined using principal component analisys methods as a
linear combination of Zernike modes, represent the best choice for wavefront reconstruction, giving us
the possibility of reducing fitting error with respect to Zernike modal expansion.
Finally, we tested the possibility to use auto-regressive moving average time series models to provide
short time forecasting for Zernike or Karhunen-Loeve coefficients. This scheme may be very helpful in
improving correction bandwidth and overcoming time delay errors in next generation solar Adaptive
Optics systems.
These preliminary results could be particularly important in developing and studying design and
performance for next European Solar Telescope.
Scarica
Ottimizzazione della ricostruzione del fronte d`onda per telescopi

Diapositiva 1

Ottimizzazione della ricostruzione del fronte d`onda per telescopi

Diapositiva 1

Le onde elettromagnetiche

buon mese di giugno 2010

Diapositiva 1

luce nè onde o particelle - Liceo Classico Scientifico XXV Aprile

Sviluppo in serie di Fourier
