% -*- mode: Noweb; noweb-code-mode: c-mode -*-

\section{Laser guide star adaptive optics}
\label{sec:lgsao}

In this section, the performance of LGS AO systems are computed.
The on--axis wavefront is eithe derived straight from the LGS or estimated with a LGS LMMSE.
All the systems employ geometric Shack--Hartmann wavefront sensors (SH--WFS).
For each system, the on--axis wavefront is estimated several times with a different random draw of the phase screens.

The components common to all the systems are defined first:
\begin{itemize}
\item the science source,
<<setup>>=
source src;
src.setup("K",ARCSEC(0) , 0, INFINITY,(N_SIDE_LENSLET+1)*(N_SIDE_LENSLET+1), "SRC");
source src_lgs;
src_lgs.setup("K",ARCSEC(0) , 0, 90e3,(N_SIDE_LENSLET+1)*(N_SIDE_LENSLET+1), "SRC_LGS");
<<cleanup>>=
src.cleanup();
src_lgs.cleanup();
@ \item the atmosphere,
<<setup>>=
atmosphere atm;
//atm.setup(20e-2,30,10e3,10,0);
atm.gmt_setup(15e-2,60);
/*
float altitude[] = {0, 10e3},
xi0[] = {0.5, 0.5},
wind_speed[] = {10, 10},
wind_direction[] = {0, 0};
atm.setup(20e-2,30,altitude, xi0, wind_speed, wind_direction);
*/
<<cleanup>>=
atm.cleanup();
@ \item the diameter of the telescope,
<<setup>>=
float D = 25; // telescope diameter in meter
@ leading to a lenslet size of:
<<setup>>=
float d = D/N_SIDE_LENSLET;
@ \item the pupil mask.
<<setup>>=
mask pupil_mask;
pupil_mask.setup( (N_SIDE_LENSLET+1)*(N_SIDE_LENSLET+1) );
<<cleanup>>=
pupil_mask.cleanup();
@ \item the wavefront sensor centroid container,
<<setup>>=
centroiding cog;
cog.setup(N_SIDE_LENSLET,1);
<<cleanup>>=
cog.cleanup();
@ \item the diffraction limited science imager,
<<setup>>=
imaging imager;
imager.setup(N_SIDE_LENSLET+1,1,4,1,1);
<<cleanup>>=
imager.cleanup();
@ \item the turbulence limited science imager,
<<setup>>=
imaging imager_turb;
imager_turb.setup(N_SIDE_LENSLET+1,1,4,1,1);
<<cleanup>>=
imager_turb.cleanup();
@ \item the science imager for LGSAO,
<<setup>>=
imaging imager_lgsao;
imager_lgsao.setup(N_SIDE_LENSLET+1,1,4,1,1);
<<cleanup>>=
imager_lgsao.cleanup();
@ \item the science imager for LGSAO LMMSE,
<<setup>>=
imaging imager_lgsao_lmmse;
imager_lgsao_lmmse.setup(N_SIDE_LENSLET+1,1,4,1,1);
<<cleanup>>=
imager_lgsao_lmmse.cleanup();
@ \item the statistical tool.
<<setup>>=
stats S;
S.setup();
<<cleanup>>=
S.cleanup();
@  \end{itemize}

@ 
The Fried geometry for a circular pupil with the intensity [[threshold]] is enforced::
<<setup>>=
float threshold = 0.5;
cog.MASK_SET = fried_geometry_setup(cog.lenslet_mask, pupil_mask.m, 
				    N_SIDE_LENSLET, 16, threshold);
@ 
The filtering properties associated with the pupil are set with:
<<setup>>=
pupil_mask.set_filter();
@ 
The science is propagated through the [[pupil_mask]] to the focal plane of the [[imager]]:
<<setup>>=
src.wavefront.M = &pupil_mask;
src.wavefront.masked();
imager.propagate(&src);
//imager.frame2file("refFrame.bin");
char plotly_name[64], plotly_dir[64];
sprintf(plotly_dir,"LGSAO/D=%.1fm, N=%d, z_lgs=%2.0fkm (%d samples)/",
        D,N_SIDE_LENSLET,src_lgs.height*1e-3,N_SAMPLE);
if (PLOTLY_LIM<N_SAMPLE) {
  SET_PLOTLY_NAME(plotly_name,plotly_dir,"frames/diffraction limited");
  imager.show_frame(plotly_name);
 }

@ A few useful variables are defined here:
<<setup>>=
int NP, NP2, k_SAMPLE;
float tau=0.;
NP = N_SIDE_LENSLET+1;
NP2 = NP*NP;
k_SAMPLE = 0;
@
The science wavefront is propagated through the atmosphere from [[src]] to [[pupil_mask]].
<<science wavefront>>=
atm.get_phase_screen(&src,d,NP,d,NP,tau);
src.wavefront.masked();
if (k_SAMPLE==PLOTLY_LIM) {
  SET_PLOTLY_NAME(plotly_name,plotly_dir,"phases/science phase screen");
  src.wavefront.show_phase(plotly_name);
 }
<<setup>>=
float science_wf_rms=0.;
<<science wavefront>>=
science_wf_rms += S.var(src.wavefront.phase, &pupil_mask, NP2);

@  and then propagated to the focal plane of the imager:
<<science wavefront>>=
imager_turb.propagate(&src);
if (k_SAMPLE==PLOTLY_LIM) {
  SET_PLOTLY_NAME(plotly_name,plotly_dir,"frames/turbulence limited");
  imager_turb.show_frame(plotly_name);
 }
@ The turbulence wavefront is saved apart for later use:
<<setup>>=
complex_amplitude phase_screen;
<<science wavefront>>=
phase_screen.setup((N_SIDE_LENSLET+1)*(N_SIDE_LENSLET+1));
phase_screen.add_phase(1,src.wavefront.phase);

@ 
All CEO programs must include the following headers which also contains the headers for all the CEO library modules.
<<main header>>=
#ifndef __CEO_H__
#include "ceo.h"
#endif
@ The size of the lenslet array is defined in the header:
<<main header>>=
#define N_SIDE_LENSLET 60
#define _N_LENSLET_ (N_SIDE_LENSLET*N_SIDE_LENSLET)
@ The number of atmosphere sample
<<main header>>=
#define N_SAMPLE 1000
#define PLOTLY_LIM (N_SAMPLE-1) // (N_SAMPLE-1) for plotting, higher for disabling plotting
@
The main function is:
<<lgsao.bin>>=
<<main header>>
void SET_PLOTLY_NAME(char *name,char *dir,char *path)
{
strcpy(name, dir);
strcat(name,path);
}
int main(int argc,char *argv[]) {
cudaSetDevice(1);
<<setup>>
fprintf(stdout," [    ] of %d samples",N_SAMPLE);
for (k_SAMPLE=0;k_SAMPLE<N_SAMPLE;k_SAMPLE++) {
  fprintf(stdout,"\r [%4d]",k_SAMPLE);
  fflush(stdout);
  atm.reset();
  <<science wavefront>>
  <<lgsao>>
  phase_screen.cleanup();
 }
 fprintf(stdout,"\n");
<<results>>
<<cleanup>>
}

@ 
The laser guide star (LGS) is given by:
<<setup>>=
source lgs;
lgs.setup("R",ARCSEC(0) , 0, src_lgs.height, (N_SIDE_LENSLET+1)*(N_SIDE_LENSLET+1), "LGS");
<<cleanup>>=
lgs.cleanup();

@  The LGS turbulence wavefront and wavefront gradient are computed with:
<<lgsao>>=
atm.get_phase_screen(&lgs,d,NP,d,NP,0);
lgs.wavefront.M = &pupil_mask;
lgs.wavefront.masked();
if (k_SAMPLE==PLOTLY_LIM) {
  SET_PLOTLY_NAME(plotly_name,plotly_dir,"phases/LGS phase screen");
  lgs.wavefront.show_phase(plotly_name);
 }
<<setup>>=
float lgs_wf_rms=0.;
<<lgsao>>=
lgs_wf_rms += S.var(lgs.wavefront.phase, &pupil_mask, NP2);
@ and with:
<<lgsao>>=
atm.get_phase_screen_gradient(&cog,N_SIDE_LENSLET,d,&lgs,tau);

@  
The LGS wavefront is subtracted from the science wavefront:
<<setup>>=
LMMSE lgs_lgs_lmmse;
lgs_lgs_lmmse.setup(&atm,&lgs,&src_lgs,d,N_SIDE_LENSLET,&pupil_mask,"MINRES");
<<lgsao>>=
lgs_lgs_lmmse.estimation(&cog);
src.wavefront.reset(phase_screen);
src.wavefront.add_phase(-1,lgs_lgs_lmmse.d__phase_est);
if (k_SAMPLE==PLOTLY_LIM) {
  SET_PLOTLY_NAME(plotly_name,plotly_dir,"phases/LGS AO");
  src.wavefront.show_phase(plotly_name);
 }
<<setup>>=
float lgs_wfe_rms=0.;
<<lgsao>>=
lgs_wfe_rms += S.var(src.wavefront.phase, 
		   &pupil_mask, NP2);
@  and the residual wavefront corresponding image is computed.
<<lgsao>>=
imager_lgsao.propagate(&src);
if (k_SAMPLE==PLOTLY_LIM) {
  SET_PLOTLY_NAME(plotly_name,plotly_dir,"frames/LGS AO");
  imager_lgsao.show_frame(plotly_name,&imager);
 }
@ 
The wavefront is reconstructed from the LGS centroids:
<<setup>>=
LMMSE lgs_lmmse;
lgs_lmmse.setup(&atm,&lgs,&src,d,N_SIDE_LENSLET,&pupil_mask,"MINRES");
int cvgce_iteration=0;
float elapsed_time=0.;
<<lgsao>>=
lgs_lmmse.estimation(&cog);
elapsed_time += lgs_lmmse.elapsed_time;
cvgce_iteration += lgs_lmmse.iSolve.cvgce_iteration;
@ The LMMMSE LGS wavefront estimate is subtracted from the science wavefront:
<<lgsao>>=
src.wavefront.reset(phase_screen);
src.wavefront.add_phase(-1,lgs_lmmse.d__phase_est);
if (k_SAMPLE==PLOTLY_LIM) {
  SET_PLOTLY_NAME(plotly_name,plotly_dir,"phases/LMMSE LGS AO");
  src.wavefront.show_phase(plotly_name);
 }
<<setup>>=
float lgs_wfe_rms_lmmse=0.;
<<lgsao>>=
lgs_wfe_rms_lmmse += S.var(src.wavefront.phase, 
		   &pupil_mask, NP2);
//lgs_lmmse.toFile("phaseEstLgsItp.bin");
<<cleanup>>=
lgs_lmmse.cleanup();
@ and the residual wavefront corresponding image is computed.
<<lgsao>>=
imager_lgsao_lmmse.propagate(&src);
//imager.frame2file("lgsLmmseFrame.bin");
if (k_SAMPLE==PLOTLY_LIM) {
  SET_PLOTLY_NAME(plotly_name,plotly_dir,"frames/LMMSE LGS AO");
  imager_lgsao_lmmse.show_frame(plotly_name,&imager);
 }

@ 
\subsection{Results}
\label{sec:results}

<<results>>=
science_wf_rms    = 1E9*sqrtf(science_wf_rms/N_SAMPLE);
lgs_wf_rms        = 1E9*sqrtf(lgs_wf_rms/N_SAMPLE);
lgs_wfe_rms       = 1E9*sqrtf(lgs_wfe_rms/N_SAMPLE);
lgs_wfe_rms_lmmse = 1E9*sqrtf(lgs_wfe_rms_lmmse/N_SAMPLE);
elapsed_time /= N_SAMPLE;
cvgce_iteration /= N_SAMPLE;
printf("------------------------------\n");
printf("___  TURBULENCE WAVEFRONT ___\n");
printf("\n NGS WF RMS: %7.2fnm\n", science_wf_rms);
printf("\n LGS WF RMS: %7.2fnm\n", lgs_wf_rms);
printf("\n___ ON-AXIS WAVEFRONT ESTIMATE FROM LGS ___\n");
printf("\n NGS-LGS WFE RMS : %8.3fnm\n",lgs_wfe_rms);
printf("\n LGS MMSE WFE RMS: %8.3fnm in %.2fms with %d iterations\n",
       lgs_wfe_rms_lmmse,elapsed_time,cvgce_iteration);
printf("------------------------------\n");