> <> int main(int argc,char *argv[]) { printf(" ##### TEST 1 #####\n"); if (argc != 4) { printf("Usage: %s N_SIDE_LENSLET N_PX_LENSLET N_GS\n",argv[0]); return 1; } cudaSetDevice(0); <> wfs.calibrate(&gs,0.5); atm.get_phase_screen(&gs,N_GS,p,N_PX,p,N_PX,0); gs.wavefront.masked(); wfs.analyze(&gs); gs.wavefront.show_amplitude("Shack-Hartmann/pupil"); gs.wavefront.show_phase("Shack-Hartmann/wavefront"); wfs.camera.show_frame("Shack-Hartmann/frame"); wfs.data_proc.show_centroids("Shack-Hartmann/centroids"); wfs.data_proc.show_flux("Shack-Hartmann/flux"); lmmse.estimation(&wfs.data_proc); atm.get_phase_screen(&gs_dm,N_GS,d,N_ACTUATOR,d,N_ACTUATOR,0); gs_dm.wavefront.masked(&DM_mask); gs_dm.wavefront.show_phase("Shack-Hartmann/wavefront low res."); phase_screen.M = &DM_mask; phase_screen.add_phase(1,lmmse.d__phase_est); //phase_screen.add_phase(-1,gs_dm.wavefront.phase); phase_screen.show_phase("Shack-Hartmann/recon. wavefront"); /* atm.get_phase_screen_gradient(&wfs.data_proc,N_SIDE_LENSLET,d,&gs,N_GS,0); wfs.data_proc.show_centroids("Shack-Hartmann/geometric centroids"); lmmse.estimation(&wfs.data_proc); phase_screen.reset(); phase_screen.add_phase(1,lmmse.d__phase_est); //phase_screen.add_phase(-1,gs_dm.wavefront.phase); phase_screen.show_phase("Shack-Hartmann/geom. recon. wavefront"); */ <> } @ The input parameters are \begin{itemize} \item the size of the lenslet array [[N_SIDE_LENSLET]]$\times$[[N_SIDE_LENSLET]], <>= int N_SIDE_LENSLET, N_ACTUATOR, N_ACTUATOR2; N_SIDE_LENSLET = atoi( argv[1] ); N_ACTUATOR = N_SIDE_LENSLET + 1; N_ACTUATOR2 = N_ACTUATOR*N_ACTUATOR; printf("__ SHACK-HARTMANN __\n"); printf(" . %dX%d lenslets: \n", N_SIDE_LENSLET, N_SIDE_LENSLET); @ \item the number of pixel per lenslet [[N_PX_LENSLET]]$\times$[[N_PX_LENSLET]] <>= int N_PX_LENSLET; N_PX_LENSLET = atoi( argv[2] ); printf(" . %dX%d pixels per lenslet: \n", N_PX_LENSLET, N_PX_LENSLET); @ \item the number of sources, <>= int N_GS; N_GS = atoi( argv[3] ); printf(" . source #: %d\n", N_GS); @ \end{itemize} The total number of pixel is given by <>= int N_PX, N_PX2; N_PX = N_SIDE_LENSLET*N_PX_LENSLET; N_PX2 = N_PX*N_PX; float p, d; p = 25.0/N_PX; d = 25.0/N_SIDE_LENSLET; @ The guide star is defined with <>= float zenith[] = {ARCSEC(30),ARCSEC(30),ARCSEC(30),ARCSEC(30),ARCSEC(30),ARCSEC(30)}, azimuth[] = {0,2.*PI/6.,4.*PI/6.,6.*PI/6.,8.*PI/6.,10.*PI/6.,12.*PI/6.}; source gs; gs.setup("K", zenith, azimuth, INFINITY, N_GS, N_PX2); <>= gs.cleanup(); @ The telescope pupil is defined with the [[mask]] structure: <>= mask telescope; telescope.setup_circular(N_PX); gs.wavefront.masked(&telescope); <>= telescope.cleanup(); @ The wavefront sampled on the actuators location is from the following guide star: <>= source gs_dm; gs_dm.setup("K", zenith, azimuth, INFINITY, N_GS, N_ACTUATOR2); <>= gs_dm.cleanup(); @ The wavefront sensor is defined with: <>= shackHartmann wfs; wfs.setup(N_SIDE_LENSLET,N_PX_LENSLET,d,2,N_PX_LENSLET,1,N_GS); <>= wfs.cleanup(); @ The timing is done with the [[stopwatch]] structure: <>= stopwatch tid; @ The square pupil is made of tiles [[N_PX_LENSLET]]$\times$[[N_PX_LENSLET]] pixels with intensity equal to lenslet index <>= dim3 blockDim(16,16); dim3 gridDim(1+N_PX_LENSLET*N_SIDE_LENSLET/16,1+N_PX_LENSLET*N_SIDE_LENSLET/16,N_GS); squareTiledPupil LLL gridDim,blockDim RRR (gs.wavefront.amplitude, N_PX_LENSLET, N_SIDE_LENSLET, N_GS); gs.wavefront.show_amplitude("Shack-Hartmann/amplitude",N_PX,N_PX*N_GS); <>= __global__ void squareTiledPupil(float* amplitude, const int N_PX_LENSLET, const int N_SIDE_LENSLET, const int N_GS) { int i, j, iLenslet, jLenslet, ij_PUPIL, iSource, kLenslet; iSource = blockIdx.z; if ( threads2lenslet(threadIdx, blockIdx, blockDim, &i, &j, N_PX_LENSLET, &iLenslet, &jLenslet, N_SIDE_LENSLET) ) { ij_PUPIL = lenslet2array(i,j,N_PX_LENSLET,iLenslet,jLenslet,N_SIDE_LENSLET,iSource); kLenslet = iLenslet*N_SIDE_LENSLET + jLenslet; amplitude[ij_PUPIL] = ((float) (kLenslet + 1)) + iSource*0.5; } } @ The atmosphere model is created with: <>= atmosphere atm; atm.gmt_setup(20e-2,30); <>= atm.cleanup(); @ The reconstructed wavefront is saved in: <>= complex_amplitude phase_screen; phase_screen.setup((N_SIDE_LENSLET+1)*(N_SIDE_LENSLET+1)); <>= phase_screen.cleanup(); @ The DM mask: <>= mask DM_mask; DM_mask.setup( (N_SIDE_LENSLET+1)*(N_SIDE_LENSLET+1) ); <>= DM_mask.cleanup(); @ The Fried geometry for a circular pupil with the intensity [[threshold]] is enforced:: <>= float threshold = 0.5; wfs.data_proc.MASK_SET = fried_geometry_setup(wfs.data_proc.lenslet_mask, DM_mask.m, N_SIDE_LENSLET, 16, threshold); @ The filtering properties associated with the pupil are set with: <>= DM_mask.set_filter(); @ The wavefront is reconstructed from the centroids: <>= LMMSE lmmse; lmmse.setup(&atm,&gs,&gs,d,N_SIDE_LENSLET,&DM_mask,"MINRES"); <>= lmmse.cleanup(); @ \subsection{Test 2} \label{sec:test-2} The input parameters are \begin{itemize} \item the size of the lenslet array [[N_SIDE_LENSLET]]$\times$[[N_SIDE_LENSLET]], <>= int N_SIDE_LENSLET, N_ACTUATOR, N_ACTUATOR2; N_SIDE_LENSLET = atoi( argv[1] ); int _N_LENSLET_ = (N_SIDE_LENSLET*N_SIDE_LENSLET); N_ACTUATOR = N_SIDE_LENSLET + 1; N_ACTUATOR2 = N_ACTUATOR*N_ACTUATOR; printf("__ SHACK-HARTMANN __\n"); printf(" . %dX%d lenslets: \n", N_SIDE_LENSLET, N_SIDE_LENSLET); @ \item the number of pixel per lenslet [[N_PX_LENSLET]]$\times$[[N_PX_LENSLET]] <>= int N_PX_LENSLET; N_PX_LENSLET = atoi( argv[2] ); printf(" . %dX%d pixels per lenslet: \n", N_PX_LENSLET, N_PX_LENSLET); @ \item the number of sources, <>= int N_GS; N_GS = atoi( argv[3] ); printf(" . source #: %d\n", N_GS); @ \end{itemize} The total number of pixel is given by <>= int N_PX, N_PX2; N_PX = N_SIDE_LENSLET*N_PX_LENSLET; N_PX2 = N_PX*N_PX; float p, d, D; D = 25.0; p = D/N_PX; d = D/N_SIDE_LENSLET; @ The components common to all the systems are defined first: \begin{itemize} \item the science source, <>= source src; src.setup("J",ARCSEC(0) , 0, INFINITY,(N_SIDE_LENSLET+1)*(N_SIDE_LENSLET+1), "SRC"); <>= src.cleanup(); @ \item the atmosphere, <>= 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); */ <>= atm.cleanup(); @ \item the pupil mask. <>= //mask pupil_mask; //pupil_mask.setup( (N_SIDE_LENSLET+1)*(N_SIDE_LENSLET+1) ); <>= //pupil_mask.cleanup(); @ \item the wavefront sensor centroid container, <>= //centroiding cog; //cog.setup(N_SIDE_LENSLET,1); <>= //cog.cleanup(); @ \item the diffraction limited science imager, <>= imaging imager; imager.setup(N_SIDE_LENSLET+1,1,4,1,1); <>= imager.cleanup(); @ \item the turbulence limited science imager, <>= imaging imager_turb; imager_turb.setup(N_SIDE_LENSLET+1,1,4,1,1); <>= imager_turb.cleanup(); @ \item the science imager for NGSAO LMMSE, <>= imaging imager_ngsao_lmmse; imager_ngsao_lmmse.setup(N_SIDE_LENSLET+1,1,4,1,1); <>= imager_ngsao_lmmse.cleanup(); @ \item the statistical tool. <>= stats S; S.setup(); <>= S.cleanup(); @ \end{itemize} The natural guide star (NGS) is given by: <>= source ngs; ngs.setup("R",ARCSEC(0) , 0, INFINITY, (N_SIDE_LENSLET+1)*(N_SIDE_LENSLET+1), "NGS"); <>= ngs.cleanup(); @ The wavefront sensor is defined with: <>= shackHartmann wfs; wfs.setup(N_SIDE_LENSLET,N_PX_LENSLET,d,2,N_PX_LENSLET,1,N_GS); wfs.calibrate(&ngs,0.5); <>= wfs.cleanup(); @ The Fried geometry for a circular pupil with the intensity [[threshold]] is enforced:: <>= //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: <>= wfs.valid_actuator.set_filter(); //pupil_mask.set_filter(); @ The science is propagated through the [[pupil_mask]] to the focal plane of the [[imager]]: <>= src.wavefront.masked(&wfs.valid_actuator); imager.propagate(&src); //imager.frame2file("refFrame.bin"); char plotly_name[64], plotly_dir[64]; sprintf(plotly_dir,"Shack-Hartmann+LMMSE/D=%.1fm, N=%d (%d samples)/",D,N_SIDE_LENSLET,N_SAMPLE); SET_PLOTLY_NAME(plotly_name,plotly_dir,"frames/diffraction limited"); imager.show_frame(plotly_name); @ A few useful variables are defined here: <>= 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]]. <>= atm.get_phase_screen(&src,d,NP,d,NP,tau); src.wavefront.masked(); if (k_SAMPLE==(N_SAMPLE-1)) { SET_PLOTLY_NAME(plotly_name,plotly_dir,"phases/science phase screen"); src.wavefront.show_phase(plotly_name); } <>= float science_wf_rms=0.; <>= science_wf_rms += S.var(src.wavefront.phase, &wfs.valid_actuator, NP2); @ and then propagated to the focal plane of the imager: <>= imager_turb.propagate(&src); if (k_SAMPLE==(N_SAMPLE-1)) { 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: <>= complex_amplitude phase_screen; 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. <

>= #ifndef __CEO_H__ #include "h" #endif @ The number of atmosphere sample <

>= #define N_SAMPLE 1 #define PLOTLY_LIM (N_SAMPLE-1) // (N_SAMPLE-1) for plotting, higher for disabling plotting @ The main function is: <>= <

> void SET_PLOTLY_NAME(char *name,char *dir,char *path) { strcpy(name, dir); strcat(name,path); } int main(int argc,char *argv[]) { printf(" ##### TEST 2 #####\n"); if (argc != 4) { printf("Usage: %s N_SIDE_LENSLET N_PX_LENSLET N_GS\n",argv[0]); return 1; } cudaSetDevice(0); <> INFO(" Samples: %d:\n",N_SAMPLE); for (k_SAMPLE=0;k_SAMPLE> <> } INFO("\n"); <> <> } @ The NGS wavefront gradients are computed with: <>= atm.get_phase_screen_gradient(&wfs.data_proc,N_SIDE_LENSLET,d,&ngs,tau); @ The wavefront is reconstructed from the NGS centroids: <>= LMMSE ngs_lmmse; ngs_lmmse.setup(&atm,&ngs,&src,&wfs,"MINRES"); int cvgce_iteration=0; float elapsed_time=0.; <>= ngs_lmmse.estimation(&wfs.data_proc); elapsed_time += ngs_lmmse.elapsed_time; cvgce_iteration += ngs_lmmse.iSolve.cvgce_iteration; @ The LMMMSE NGS wavefront estimate is subtracted from the science wavefront: <>= src.wavefront.reset(phase_screen); src.wavefront.add_phase(-1,ngs_lmmse.d__phase_est); if (k_SAMPLE==PLOTLY_LIM) { SET_PLOTLY_NAME(plotly_name,plotly_dir,"phases/LMMSE NGS AO"); src.wavefront.show_phase(plotly_name); } <>= float ngs_wfe_rms_lmmse=0.; <>= ngs_wfe_rms_lmmse += S.var(src.wavefront.phase, &wfs.valid_actuator, NP2); //ngs_lmmse.toFile("phaseEstNgsItp.bin"); <>= ngs_lmmse.cleanup(); @ and the residual wavefront corresponding image is computed. <>= imager_ngsao_lmmse.propagate(&src); //imager.frame2file("ngsLmmseFrame.bin"); if (k_SAMPLE==PLOTLY_LIM) { SET_PLOTLY_NAME(plotly_name,plotly_dir,"frames/LMMSE NGS AO"); imager_ngsao_lmmse.show_frame(plotly_name,&imager); } @ \subsubsection{Results} \label{sec:results} <>= science_wf_rms = 1E9*sqrtf(science_wf_rms/N_SAMPLE); ngs_wfe_rms_lmmse = 1E9*sqrtf(ngs_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___ ON-AXIS WAVEFRONT ESTIMATE FROM NGS ___\n"); printf("\n NGS MMSE WFE RMS: %8.3fnm in %.2fms with %d iterations\n", ngs_wfe_rms_lmmse,elapsed_time,cvgce_iteration); printf("------------------------------\n"); @ \subsection{Python} \label{sec:python-1} <>= #!/usr/bin/env python import sys sys.path.append("/home/rconan/CEO/python") from pylab import * from import * nLenslet = 20 D = 10.0 n = 7 nPx = n*nLenslet nn = array([nPx,nPx],dtype=int32) nSrc = 1 zen = zeros( nSrc, dtype=float32) azi = zeros( nSrc, dtype=float32) src = Source("K",zen,azi,float("inf"),nn,fwhm=2) tel = Telescope(n*nLenslet) src.masked(tel) wfs = ShackHartmann(nLenslet, n, D/nLenslet, 2, n, 1, nSrc) wfs.calibrate(src,0.5) wfs.analyze(src) frame0 = zeros( (nPx*nSrc,nPx), order="c", dtype=float32) wfs.getframe(frame0) figure(1) ax = imshow(frame0,interpolation='None') colorbar(ax) draw() atm = GmtAtmosphere(15e-2,30) p = D/nPx atm.get_phase_screen(src, nSrc, p, nPx, p, nPx, 0.0) src.masked(tel) ps = zeros( (nPx*nSrc,nPx), order="c", dtype=float32) src.getphase(ps) figure(2) subplot(1, 3, 1) ax = imshow(ps,interpolation='None') colorbar(ax) wfs.reset() wfs.analyze(src) wfs.getframe(frame0) subplot(1, 3, 2) ax = imshow(frame0,interpolation='None') colorbar(ax) c = zeros( nSrc*2*nLenslet**2, order="c", dtype=float32 ) wfs.getc(c) subplot(1, 3, 3) ax = imshow(reshape(c,[nSrc*nLenslet*2,nLenslet]),interpolation='None') colorbar(ax) show()