% -*- mode: Noweb; noweb-code-mode: c-mode -*- @ \index{gmtMirrors} \section{The files} \label{sec:files} \subsection{Header} \label{sec:header} <>= #ifndef __GMTMIRRORS_H_ #define __GMTMIRRORS_H_ #ifndef __RAYTRACING_H__ #include "rayTracing.h" #endif #ifndef __SOURCE_H__ #include "source.h" #endif __device__ rtd bending_modes_surface(vector *v, const int N_mode, rtd *b, double *BM, const int BM_N_SAMPLE, const double BM_radius); __device__ rtd partial_x_bending_modes_surface(vector *v, const int N_mode, rtd *b, double *BM, const int BM_N_SAMPLE, const double BM_radius); __device__ rtd partial_y_bending_modes_surface(vector *v, const int N_mode, rtd *b, double *BM, const int BM_N_SAMPLE, const double BM_radius); __device__ rtd partial_z_bending_modes_surface(void); <> <> <> <> <> <> #endif // __GMTMIRRORS_H_ @ \subsection{Source} \label{sec:source} <>= #include "gmtMirrors.h" <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> @ @ \section{GMT M1} \label{sec:gmt-m1} \index{gmtMirrors!gmt\_m1} A new structure to hold GMT M1 and M2 parameters and functions is defined. <>= struct gmt_m1 { <> <> <> void preset(bundle *rays, rtd margin); void edge_sensors_data(void); }; @ with <>= void setup(void); //void setup(modes *_BS_); void setup(char *_filename_, int _N_, int _n_mode_); void setup(zernikeS *ZS); void cleanup(void); void update(vector _origin_, vector _euler_angles_,int idx); void reset(void); void trace(bundle *rays); void traceall(bundle *rays); void blocking(bundle *rays); void global_tiptilt(float tip, float tilt); void test_ray_tracing(void); void track(float *x, float *y, float *z, int N, int idx); void locate(float *x, float *y, float *z, int N, int idx); void remove(int *seg_ID, int N_ID); void keep(int *seg_ID, int N_ID); void update_conic_c(rtd *_conic_c_); void update_conic_k(rtd *_conic_k_); void set_reflectivity(float *reflectivity); @ The GMT M12 parameters are: \begin{itemize} \item the mirror ID: <>= int M_ID; @ \item the mirror assembly diameter $[[D_assembly]]$: <>= rtd D_assembly; @ \item the segment clear aperture diameter $[[D_clear]]$: <>= rtd D_clear; @ \item the segment full aperture diameter $[[D_full]]$: <>= rtd D_full; @ \item the obscuration ratio of the center segment $[[ri]]$ <>= rtd ri; @ \item the tilt angle of the peripheral segment $[[beta]]$degree <>= rtd beta; @ \item the distance from the optical axis to the center of the tilted peripheral segments $[[L]]$m: <>= rtd L; @ \item the size of the square array containing M1 in meter and pixel [[D]] and [[D_px]]: <>= //rtd D; //int D_px; @ \item the mirror total area as seen by an on--axis source in square meter and number of pixel <>= rtd area0, area_fraction; float area0_px; @ \item the mirror total area as seen by the star propagating through the telescope <>= rtd area; @ \item the number of segment $[[N]]=7$ <>= int N; @ \item the edge to center depth of the segment conic $[[depth]]$: <>= rtd depth; @ \item the reference frames of the apertures: <>= coordinate_system aperture_CS; @ \item the reference frames and origins of the conics: <>= coordinate_system conic_CS; vector conic_origin[7], *d__conic_origin; @ \item the conic parameters: <>= rtd conic_c, conic_k; rtd *d__conic_c, *d__conic_k; @ \item the reference frames of the rigid bodies: <>= coordinate_system rigid_body_CS; @ \item the reference frames of the motions: <>= coordinate_system motion_CS; @ \item the mirror height <>= rtd height; @ \item the M1 pupil mask [[V]]: <>= mask *V; @ \item the segment index offset (0 for M1 and 3 for M2): <>= int idx_offset; @ \item a Zernike surface [[ZS]]: <>= zernikeS *ZS; @ \item the segment piston mask <>= int *d__piston_mask; @ \item global tip--tilt transformation variables: \begin{itemize} \item a coordinate system: <>= coordinate_system TT_CS; @ \item a pointer to a matrix array: <>= double *d__C, *d__D; @ \item a CUBLAS handle: <>= cublasHandle_t handle; @ \end{itemize} \item a 7 element array where 1 means the segment is active and 0 means the segment is missing <>= char *d__valid_segments; @ \item a bending modes structure [[BS]]: <>= //bending_modes *BS; modes BS; @ \item a 7 element array for the segment reflectivity coefficients in the range $[0,1]$ <>= float *d__segment_reflectivity; @ \end{itemize} \subsection{Setup \& Cleanup} \label{sec:setup--cleanup-1} The GMT M1 segment have all the same clear aperture of 8.365m diameter. The center segment has a hole of 2.4412m diameter. The center aperture is centered on the vertex of the conic surface that defines M1. The vertex is set as the origin of the global coordinate system (GCS). The center of the peripheral segments are evenly located on a circle of radius 8.710m with a phase of 30 degrees. The peripheral segments also tilted inwards by 13.522 degrees. \index{gmtMirrors!gmt\_m1!setup} <>= void gmt_m1::setup(void) { // BS = NULL; ZS = NULL; <> <> } @ <>= void gmt_m1::setup(modes *_BS_) { BS = _BS_; ZS = NULL; <> <> } @ <>= void gmt_m1::setup(char *_filename_, int _N_, int _n_mode_) { BS.setup(_filename_, _N_, _n_mode_); ZS = NULL; <> <> } @ or \index{gmtMirrors!gmt\_m1!setup} <>= void gmt_m1::setup( zernikeS *_ZS_) { ZS = _ZS_; //BS = NULL; <> <> } @ where <>= M_ID = 1; D_assembly = 25.498; D_full = 8.417; D_clear = 8.365; area = 357; area0 = 357; ri = 2.875/8.417; beta = 13.601685*PI/180.0; L = 8.710; N = 7; conic_c = 1.0/36.0; conic_k = 1-0.9982857; height = 0.0; idx_offset = 0; V = NULL; @ and <>= int n_byte; float segment_reflectivity[7] = {1.,1.,1.,1.,1.,1.,1.}; n_byte = 7*sizeof(float); HANDLE_ERROR( cudaMalloc((void**)&d__segment_reflectivity, n_byte ) ); HANDLE_ERROR( cudaMemcpy( d__segment_reflectivity, segment_reflectivity, n_byte, cudaMemcpyHostToDevice ) ); rtd *h__conic_c; n_byte = 7*sizeof(rtd); h__conic_c = (rtd *)malloc( n_byte ); for (int k=0; k<7; k++) h__conic_c[k] = conic_c; HANDLE_ERROR( cudaMalloc((void**)&d__conic_c, n_byte ) ); HANDLE_ERROR( cudaMemcpy( d__conic_c, h__conic_c, n_byte, cudaMemcpyHostToDevice ) ); free(h__conic_c); rtd *h__conic_k; n_byte = 7*sizeof(rtd); h__conic_k = (rtd *)malloc( n_byte ); for (int k=0; k<7; k++) h__conic_k[k] = conic_k; HANDLE_ERROR( cudaMalloc((void**)&d__conic_k, n_byte ) ); HANDLE_ERROR( cudaMemcpy( d__conic_k, h__conic_k, n_byte, cudaMemcpyHostToDevice ) ); free(h__conic_k); area_fraction = 1.0; vector __v0(0.0,0.0,0.0), __v(D_full*0.5,0.0,0.0), origin[12], euler_angles[12]; rtd D_c, o, zo; int k, idx; char tag[16]; depth = conic_equation(&__v,&__v0,conic_k,conic_c, 0, &D_c); D_c = L; __v.x = L; zo = conic_equation(&__v,&__v0,conic_k,conic_c, 0, &D_c); HANDLE_ERROR( cudaMalloc((void**)&d__valid_segments, sizeof(char)*7) ); HANDLE_ERROR( cudaMemset( d__valid_segments, 1, sizeof(char)*7) ); <> <> <> <> <> @ \subsubsection{Rigid body coordinate systems} \label{sec:rigid-body-coord} The rigid body coordinate systems are aligned with the segment coordinate systems. The origins $[x_O, y_O, z_O]$ of the segment coordinate systems are given with respect to the global coordinate system which origin is at the vertex of M1 conic surface: \begin{eqnarray} \label{eq:1} x_{O,k} &=& L\cos(\theta_{O,k}), \\ y_{O,k} &=& L\sin(\theta_{O,k}), \\ z_O &=& F(L),\\ \end{eqnarray} with, for M1, $L = 8.71$m, $\theta_{O,k}=\pi(3-2k)/6, \forall k \in [0,5]$, $x_{O,6}=y_{O,6}=0$ and \begin{equation} \label{eq:2} F(\rho) = {c\rho^2 \over 1 + \sqrt{1 - \kappa c^2 \rho^2} }. \end{equation} In addition to the change of origin with respect to the GCS, the segment coordinate system is rotated around the z--axis of $-k\pi/3$radian and then tilted around the x--axis of $\beta$. \begin{center} \input{rigidBodyCS.tex} \end{center} <>= origin[N-1].x = origin[N-1].y = 0.0; origin[N-1].z = height; euler_angles[N-1].x = 0.0; euler_angles[N-1].y = (M_ID==1) ? 0.0 : PI; euler_angles[N-1].z = (idx_offset/2)*PI; for (k=0; k> origin[idx].x = D_c*cos(o); origin[idx].y = D_c*sin(o); origin[idx].z = height + zo; if (M_ID==1) { euler_angles[idx].x = beta; euler_angles[idx].y = 0.0; euler_angles[idx].z = -PI*k/3.0; } if (M_ID==2) { euler_angles[idx].x = -beta; euler_angles[idx].y = PI; euler_angles[idx].z = -PI*k/3.0; } } sprintf(tag,"M%d RIGID BODY",M_ID); rigid_body_CS.setup(origin, euler_angles, N, tag); @ with <>= idx = (k + idx_offset)%6; @ \subsubsection{Aperture coordinate systems} \label{sec:apert-coord-syst} The coordinate systems of the segment apertures are defined with respect to the coordinate systems of the rigid body motions. The apertures are directly above the segments such as they rest on the rim of the segments. <>= origin[N-1].x = origin[N-1].y = 0.0; origin[N-1].z = (M_ID==1) ? depth: -depth; euler_angles[N-1].x = euler_angles[N-1].y = euler_angles[N-1].z = 0.0; for (k=0; k> origin[idx].x = 0.0; origin[idx].y = 0.0; origin[idx].z = (M_ID==1) ? depth: -depth; euler_angles[idx].x = 0.0; euler_angles[idx].y = 0.0; euler_angles[idx].z = 0.0; } sprintf(tag,"M%d APERTURE",M_ID); aperture_CS.setup(origin, euler_angles, N, tag); @ \subsubsection{Conic coordinate systems} \label{sec:conic-coord-syst} The coordinate systems of the segment conics are defined with respect to the coordinate systems of the rigid body motions. Both coordinate systems share the same origin but x--axis and y--axis are parallel to the x--axis and y--axis of the GCS. The origins of the segment conics are also given in the coordinate systems of the segment conics. \begin{center} \begin{tikzpicture} \draw[->] (-25mm,0) -- (35mm,0) node[below] {$\bar x$}; \draw[->] (0,-25mm) -- (0,35mm) node[right] {$\bar y$}; \draw[thin,dashed] (0,0) circle [radius=21mm]; \coordinate (O) at (0,0); \draw (O) circle [radius=10mm] node[above left=3mm] {7}; \draw[->,thick,red] ($(O)+(-5mm,0)$) -- ($(O)+(5mm,0)$) node[below] {$\hat x$}; \draw[->,thick,red] ($(O)+(0,-5mm)$) -- ($(O)+(0,5mm)$) node[right] {$\hat y$}; \foreach \x in {1,...,6} { \coordinate (O) at (150-\x*60:21mm); \draw (O) circle [radius=10mm] node[above left=3mm] {\x}; \draw[->,thick,red] ($(O)+(-5mm,0)$) -- ($(O)+(5mm,0)$) node[below] {$\check x$}; \draw[->,thick,red] ($(O)+(0,-5mm)$) -- ($(O)+(0,5mm)$) node[right] {$\check y$}; } \end{tikzpicture} \end{center} <>= origin[N-1].x = origin[N-1].y = origin[N-1].z = 0.0; euler_angles[N-1].x = 0.0; euler_angles[N-1].y = (M_ID==1) ? 0.0 : PI; euler_angles[N-1].z = (idx_offset/2)*PI; conic_origin[N-1].x = conic_origin[N-1].y = conic_origin[N-1].z = 0.0; for (k=0; k> o = PI*(3-2*k)/6.0; origin[idx].x = 0.0; origin[idx].y = 0.0; origin[idx].z = 0.0; if (M_ID==1) { euler_angles[idx].x = beta; euler_angles[idx].y = 0.0; euler_angles[idx].z = -PI*k/3.0; } if (M_ID==2) { euler_angles[idx].x = -beta; euler_angles[idx].y = PI; euler_angles[idx].z = -PI*k/3.0; } conic_origin[idx].x = -D_c*cos(o); conic_origin[idx].y = -D_c*sin(o); conic_origin[idx].z = -zo; } //euler_angles[0].x = 30.0*PI/180.0/3600.0; sprintf(tag,"M%d CONIC",M_ID); conic_CS.setup(origin, euler_angles, N, tag); HANDLE_ERROR( cudaMalloc((void**)&d__conic_origin, sizeof(vector)*N ) ); HANDLE_ERROR( cudaMemcpy( d__conic_origin, conic_origin, sizeof(vector)*N, cudaMemcpyHostToDevice ) ); @ \subsubsection{Segment motion coordinate systems} \label{sec:segm-moti-coord} The motion coordinate systems are defined with respect to the coordinate systems of the rigid body motions. For unperturbed segments, both coordinate systems coincide. <>= for (k=0; k>= void gmt_m1::cleanup(void) { INFO("@(CEO)>gmt_m1: freeing memory!\n"); <> } @ with <>= INFO(" |-"); BS.cleanup(); INFO(" |-"); <> INFO(" |-"); aperture_CS.cleanup(); INFO(" |-"); conic_CS.cleanup(); INFO(" |-"); rigid_body_CS.cleanup(); INFO(" |-"); motion_CS.cleanup(); HANDLE_ERROR( cudaFree( d__conic_origin ) ); HANDLE_ERROR( cudaFree( d__valid_segments ) ); HANDLE_ERROR( cudaFree( d__conic_c ) ); HANDLE_ERROR( cudaFree( d__conic_k ) ); HANDLE_ERROR( cudaFree( d__segment_reflectivity ) ); @ The segment radius of curvature is updated with: \index{gmtMirrors!gmt\_m1!update\_conic\_c} <>= void gmt_m1::update_conic_c(rtd *_conic_c_) { HANDLE_ERROR( cudaMemcpy( d__conic_c, _conic_c_, sizeof(rtd)*7, cudaMemcpyHostToDevice ) ); } @ The segment conic constant is updated with: \index{gmtMirrors!gmt\_m1!update\_conic\_k} <>= void gmt_m1::update_conic_k(rtd *_conic_k_) { HANDLE_ERROR( cudaMemcpy( d__conic_k, _conic_k_, sizeof(rtd)*7, cudaMemcpyHostToDevice ) ); } @ The rigid body motion parameters are updated with: \index{gmtMirrors!gmt\_m1!update} <>= void gmt_m1::update(vector _origin_, vector _euler_angles_,int idx) { <> } @ where <>= //fprintf(stdout,"\n\x1B[31m"); motion_CS.update( _origin_, _euler_angles_, --idx); //motion_CS.info(); @ The rigid body motion are reset to 0 with: \index{gmtMirrors!gmt\_m1!reset} <>= void gmt_m1::reset(void) { <> } @ where <>= vector zero; zero.x = 0.0; zero.y = 0.0; zero.z = 0.0; //fprintf(stdout,"\n\x1B[31m"); for (int k=0; k>= void gmt_m1::preset(bundle *rays, rtd margin) { dim3 blockDim(N_THREAD,N_THREAD); dim3 gridDim(rays->N_L/N_THREAD+1,rays->N_L/N_THREAD+1, rays->N_BUNDLE); m1_preset_kernel LLL gridDim , blockDim RRR (rays->d__ray, rays->N_RAY, rays->L, rays->N_L, 0.5*(margin+D_assembly)); } @ with the kernel <>= __global__ void m1_preset_kernel(ray *d__ray, int N_RAY, rtd L, int N_L, rtd R) { int i, j, k, iSource; rtd x, y, r; i = blockIdx.x * blockDim.x + threadIdx.x; j = blockIdx.y * blockDim.y + threadIdx.y; iSource = blockIdx.z; if ( (iR) d__ray[k].v = 0; } } @ \subsection{Blocking} \label{sec:blocking} The rays blocking with M1 is computed with \index{gmtMirrors!gmt\_m1!blocking} <>= void gmt_m1::blocking(bundle *rays) { <> } @ where <>= rtd R2, Rri2; R2 = D_clear*D_clear*0.25; Rri2 = R2*ri*ri; V = &(rays->V); // fprintf(stdout,"R2=%5.2f - Rri2=%5.2f\n",R2,Rri2); // printf("N_BUNDLE=%d\n",rays->N_BUNDLE); int nel = rays->N_RAY*rays->N_BUNDLE; dim3 blockDim(N_THREAD2); dim3 gridDim(nel/N_THREAD2+1); fill_ones_char LLL gridDim,blockDim RRR (V->m,nel); blockDim = dim3(N_THREAD,N_THREAD); gridDim = dim3(rays->N_RAY/N_THREAD2+1,1,rays->N_BUNDLE); m1_blocking_kernel LLL gridDim , blockDim RRR (V->m, rays->d__ray, rays->N_RAY, Rri2, R2, aperture_CS.d__R, aperture_CS.d__origin, rigid_body_CS.d__R, rigid_body_CS.d__origin, motion_CS.d__R, motion_CS.d__origin); intersection LLL gridDim , blockDim RRR (V->m, rays->d__ray, rays->N_RAY); @ and with the kernel <>= __global__ void m1_blocking_kernel(char *mask, ray *d__ray, int N_RAY, rtd inner2, rtd outer2, rtd *d__aperture_R, vector *d__aperture_origin, rtd *d__rigid_body_R, vector *d__rigid_body_origin, rtd *d__motion_R, vector *d__motion_origin) { int i, j, ij, iCoordSys, iSource; rtd rho2; rtd x, y, z, x1, y1, s0, k, l, m; vector xyz, klm, xyz_GS, klm_GS; i = blockIdx.x * blockDim.x + threadIdx.x; j = threadIdx.y; // iCoordSys = blockIdx.y; ij = j * gridDim.x * blockDim.x + i; iSource = blockIdx.z; j = ij; ij += iSource*N_RAY; if ( ( j> <> if (rho2<=outer2) { mask[ij] = 0; } } } } @ \subsection{Ray tracing} \label{sec:ray-tracing} \begin{figure} \centering \input{raytracingWorkflow.tex} \caption{Ray tracing flowchart for GMT M1 and M2 mirrors.} \label{fig:1} \end{figure} The rays propagation through M1 is computed with \index{gmtMirrors!gmt\_m1!trace} <>= void gmt_m1::trace(bundle *rays) { rays->V.area = area0*rays->N_BUNDLE*area_fraction; dim3 blockDim, gridDim; rtd R2, Rri2; R2 = D_clear*D_clear*0.25; Rri2 = R2*ri*ri; <> } @ where <>= /* if (V!=NULL) { */ /* blockDim = dim3(N_THREAD,N_THREAD); */ /* gridDim = dim3(rays->N_RAY/N_THREAD2+1,1,rays->N_BUNDLE); */ /* intersection LLL gridDim , blockDim RRR (V->m, rays->d__ray, rays->N_RAY); */ /* } */ V = &(rays->V); // fprintf(stdout,"R2=%5.2f - Rri2=%5.2f\n",R2,Rri2); blockDim = dim3(1,1); gridDim = dim3(1,1,rays->N_BUNDLE); m1_trace_chief_kernel LLL gridDim , blockDim RRR (V->m, rays->d__chief_ray, 1, conic_CS.d__R, conic_CS.d__origin, conic_k, conic_c, d__conic_origin, rigid_body_CS.d__R, rigid_body_CS.d__origin, motion_CS.d__R, motion_CS.d__origin); blockDim = dim3(N_THREAD,N_THREAD); gridDim = dim3(rays->N_RAY/N_THREAD2+1,1,rays->N_BUNDLE); HANDLE_ERROR( cudaMemset(V->m, 0, sizeof(char)*rays->N_RAY*rays->N_BUNDLE ) ); HANDLE_ERROR( cudaMemset(rays->d__piston_mask, 0, sizeof(int)*rays->N_RAY*rays->N_BUNDLE ) ); /*if (ZS!=NULL) m1_trace_kernel LLL gridDim , blockDim RRR (V->m, rays->d__ray, rays->N_RAY, Rri2, R2, aperture_CS.d__R, aperture_CS.d__origin, conic_CS.d__R, conic_CS.d__origin, d__conic_k, d__conic_c, d__conic_origin, rigid_body_CS.d__R, rigid_body_CS.d__origin, motion_CS.d__R, motion_CS.d__origin, ZS->max_n, ZS->n_mode, ZS->d__a, ZS->d__cx, ZS->d__cy, 0.5*D_full, rays->d__chief_ray, rays->d__piston_mask, d__valid_segments, d__segment_reflectivity); if (BS!=NULL)*/ m1_bm_trace_kernel LLL gridDim , blockDim RRR (V->m, rays->d__ray, rays->N_RAY, Rri2, R2, aperture_CS.d__R, aperture_CS.d__origin, conic_CS.d__R, conic_CS.d__origin, d__conic_k, d__conic_c, d__conic_origin, rigid_body_CS.d__R, rigid_body_CS.d__origin, motion_CS.d__R, motion_CS.d__origin, BS.n_mode, BS.d__b, BS.d__BMS, BS.BM_N_SAMPLE, BS.BM_radius, rays->d__chief_ray, rays->d__piston_mask, d__valid_segments, d__segment_reflectivity); intersection LLL gridDim , blockDim RRR (V->m, rays->d__ray, rays->N_RAY); float previous_nnz = rays->V.nnz; //fprintf(stdout,"M%d previous nnz rays: %f\n",M_ID,previous_nnz); rays->V.set_filter_quiet(); //fprintf(stdout,"M%d current nnz rays: %f\n",M_ID,rays->V.nnz); if (previous_nnzN_RAY_TOTAL) rays->V.area *= rays->V.nnz/previous_nnz; //fprintf(stdout,"M%d current area: %f\n",M_ID,rays->V.area); @ and with the kernel <>= __global__ void m1_trace_kernel(char *mask, ray *d__ray, int N_RAY, rtd inner2, rtd outer2, rtd *d__aperture_R, vector *d__aperture_origin, rtd *d__conic_R, vector *d__conic_origin, rtd *d__Fk, rtd *d__Fc, vector *d__conic_self_origin, rtd *d__rigid_body_R, vector *d__rigid_body_origin, rtd *d__motion_R, vector *d__motion_origin, int max_n, int n_mode, rtd *d__a, rtd *d__cx, rtd *d__cy, rtd R, ray *d__chief_ray, int *d__piston_mask, char *d__valid_segments, float *reflectivity) { int i, j, ij, iCoordSys, iSource, b_n_mode; rtd rho2; rtd x, y, z, x1, y1, s0, k, l, m; rtd s1, S, K, L ,M, dSds; rtd G2, a; rtd Fc, Fk; vector d__origin, xyz, klm, xyz_GS, klm_GS, Snormal; i = blockIdx.x * blockDim.x + threadIdx.x; j = threadIdx.y; // iCoordSys = blockIdx.y; ij = j * gridDim.x * blockDim.x + i; iSource = blockIdx.z; j = ij; ij += iSource*N_RAY; b_n_mode = max_n*(max_n+1)*0.5; if ( ( j0) { d__origin = d__conic_self_origin[iCoordSys]; <> <> if (rho2<=outer2) { mask[ij] = 1; d__piston_mask[ij] = iCoordSys + 1; <> d__ray[ij].optical_path_difference += d__ray[ij].optical_path_length - d__chief_ray[iSource].optical_path_length; d__ray[ij].throughput *= reflectivity[iCoordSys]; return; } } } } } @ for the bending modes, with the kernel <>= __global__ void m1_bm_trace_kernel(char *mask, ray *d__ray, int N_RAY, rtd inner2, rtd outer2, rtd *d__aperture_R, vector *d__aperture_origin, rtd *d__conic_R, vector *d__conic_origin, rtd *d__Fk, rtd *d__Fc, vector *d__conic_self_origin, rtd *d__rigid_body_R, vector *d__rigid_body_origin, rtd *d__motion_R, vector *d__motion_origin, const int N_mode, rtd *b, double *BM, const int BM_N_SAMPLE, const double BM_radius, ray *d__chief_ray, int *d__piston_mask, char *d__valid_segments, float *reflectivity) { int i, j, ij, iCoordSys, iSource; rtd rho2; rtd x, y, z, x1, y1, s0, k, l, m; rtd s1, S, K, L ,M, dSds; rtd G2, a; rtd Fc, Fk; vector d__origin, xyz, klm, xyz_GS, klm_GS, Snormal; i = blockIdx.x * blockDim.x + threadIdx.x; j = threadIdx.y; // iCoordSys = blockIdx.y; ij = j * gridDim.x * blockDim.x + i; iSource = blockIdx.z; j = ij; ij += iSource*N_RAY; if ( ( j0) { d__origin = d__conic_self_origin[iCoordSys]; <> <> if (rho2<=outer2) { mask[ij] = 1; d__piston_mask[ij] = iCoordSys + 1; <> d__ray[ij].optical_path_difference += d__ray[ij].optical_path_length - d__chief_ray[iSource].optical_path_length; d__ray[ij].throughput *= reflectivity[iCoordSys]; return; } } } } } @ and the chief ray kernel <>= __global__ void m1_trace_chief_kernel(char *mask, ray *d__ray, int N_RAY, rtd *d__conic_R, vector *d__conic_origin, const rtd Fk, const rtd Fc, vector *d__conic_self_origin, rtd *d__rigid_body_R, vector *d__rigid_body_origin, rtd *d__motion_R, vector *d__motion_origin) { int j, ij, iCoordSys, iSource; rtd x1, y1, s0, k, l, m; rtd s1, S, K, L ,M, dSds; rtd G2, a; vector d__origin, xyz, klm, xyz_GS, klm_GS; iSource = blockIdx.z; ij = iSource*N_RAY; iCoordSys = 0; d__origin = d__conic_self_origin[iCoordSys]; <> <> <> <> <> <> } @ The following is a test routine to validate the sequence of transformation in the ray tracing kernel: <>= void gmt_m1::test_ray_tracing(void) { int iCoordSys; vector in, out; in.x = 8.71*cos(PI/6); in.y = 8.71*sin(PI/6); in.z = 0.0; printf("Output vector:\n"); for (iCoordSys = 0; iCoordSys>= <> <> @ (i) the rays are transformed into the CS of the segment, <>= // RIGID BODY >>> forward_transform(&xyz_GS, &(d__ray[ij].coordinates), d__rigid_body_R+iCoordSys*9, d__rigid_body_origin+iCoordSys); forward_transform_centered(&klm_GS, &(d__ray[ij].directions), d__rigid_body_R+iCoordSys*9); @ (ii) the segments are perturbed <>= // MOTION >>> forward_transform(&xyz_GS, &xyz_GS, d__motion_R+iCoordSys*9, d__motion_origin+iCoordSys); forward_transform_centered(&klm_GS, &klm_GS, d__motion_R+iCoordSys*9); @ and (iii) the rays are transformed into the CS of the aperture: <>= // APERTURE >>> forward_transform(&xyz, &xyz_GS, d__aperture_R+iCoordSys*9, d__aperture_origin+iCoordSys); forward_transform_centered(&klm, &klm_GS, d__aperture_R+iCoordSys*9); x = xyz.x; y = xyz.y; z = xyz.z; k = klm.x; l = klm.y; m = klm.z; <> rho2 = x1*x1 + y1*y1; if ( (iCoordSys==6) && (rho2>= <> <> <> <> <> <> @ \begin{itemize} \item Zernike modes: <>= klm.x = klm_GS.x; klm.y = klm_GS.y; klm.z = klm_GS.z; xyz.x = xyz_GS.x; xyz.y = xyz_GS.y; xyz.z = xyz_GS.z; <> <> klm_GS.x = klm.x; klm_GS.y = klm.y; klm_GS.z = klm.z; xyz_GS.x = xyz.x; xyz_GS.y = xyz.y; xyz_GS.z = xyz.z; <> <> @ \item bending modes: <>= klm.x = klm_GS.x; klm.y = klm_GS.y; klm.z = klm_GS.z; xyz.x = xyz_GS.x; xyz.y = xyz_GS.y; xyz.z = xyz_GS.z; <> <> klm_GS.x = klm.x; klm_GS.y = klm.y; klm_GS.z = klm.z; xyz_GS.x = xyz.x; xyz_GS.y = xyz.y; xyz_GS.z = xyz.z; <> <> @ \end{itemize} <>= // CONIC >>> backward_transform_centered(&xyz, &xyz_GS, d__conic_R+iCoordSys*9); backward_transform_centered(&klm, &klm_GS, d__conic_R+iCoordSys*9); @ <>= G2 = K*K + L*L + M*M; a = k*K + l*L + m*M; a *= -2.0/G2; klm.x += a*K; klm.y += a*L; klm.z += a*M; @ <>= G2 = K*K + L*L + M*M; a = klm_GS.x*K + klm_GS.y*L + klm_GS.z*M;x a *= -2.0/G2; klm_GS.x += a*K; klm_GS.y += a*L; klm_GS.z += a*M; @ <>= // CONIC <<< forward_transform_centered(&xyz_GS, &xyz, d__conic_R+iCoordSys*9); forward_transform_centered(&klm_GS, &klm, d__conic_R+iCoordSys*9); @ <>= // MOTION <<< backward_transform(&xyz_GS, &xyz_GS, d__motion_R+iCoordSys*9, d__motion_origin+iCoordSys); backward_transform_centered(&klm_GS, &klm_GS, d__motion_R+iCoordSys*9); @ <>= // RIGID_BODY <<< backward_transform(&(d__ray[ij].coordinates), &xyz_GS, d__rigid_body_R+iCoordSys*9, d__rigid_body_origin+iCoordSys); backward_transform_centered(&(d__ray[ij].directions), &klm_GS, d__rigid_body_R+iCoordSys*9); @ The vignetting by the aperture is computed as the intersection between the aperture mask and the ray vignetting flags: <>= __global__ void intersection(char *mask, ray *d__ray, int N_RAY) { int i, j, ij, iSource; i = blockIdx.x * blockDim.x + threadIdx.x; j = threadIdx.y; ij = j * gridDim.x * blockDim.x + i; iSource = blockIdx.z; if ( ij>= if (m==0) { return; } s0 = -z/m; x1 = x + k*s0; y1 = y + l*s0; @ The intersection with the plane $z=0$ in the conic coordinate transformed is computed with \begin{itemize} \item Zernike modes: <>= k = klm.x; l = klm.y; m = klm.z; if (m==0) { return; } s0 = -xyz.z/m; x1 = xyz.x + k*s0; y1 = xyz.y + l*s0; d__ray[ij].optical_path_length = s0; s0 = s1 = 0; for (j=0; j> dSds = K*k + L*l + M*m; if (dSds==0) { break; } s1 = s0 - S/dSds; if (abs(s1-s0)> d__ray[ij].optical_path_length += s1; d__ray[ij].n_iteration = j; break; } s0 = s1; } @ \item bending modes: <>= k = klm.x; l = klm.y; m = klm.z; if (m==0) { return; } s0 = -xyz.z/m; x1 = xyz.x + k*s0; y1 = xyz.y + l*s0; d__ray[ij].optical_path_length = s0; s0 = s1 = 0; for (j=0; j> dSds = K*k + L*l + M*m; if (dSds==0) { break; } s1 = s0 - S/dSds; if (abs(s1-s0)> d__ray[ij].optical_path_length += s1; d__ray[ij].n_iteration = j; break; } s0 = s1; } @ \end{itemize} The aspheric surface is the sum of a conic and an aspheric surface (Zernike of bending modes). The conic surface is defined in the conic CS whereas the aspheric surface is defined in the motion CS. The coordinates are transformed into the conic CS where the surface height is computed: <>= // CONIC >>> backward_transform_centered(&xyz, &xyz, d__conic_R+iCoordSys*9); S = conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x)); @ The coordinates are transformed back into the motion CS where the Zernike height is added to the conic height: \begin{itemize} \item Zernike modes: <>= // CONIC <<< forward_transform_centered(&xyz, &xyz, d__conic_R+iCoordSys*9); S += zernike_surface(&xyz, &d__origin, max_n, d__a + iCoordSys*n_mode, R); @ \item bending modes: <>= // CONIC <<< forward_transform_centered(&xyz, &xyz, d__conic_R+iCoordSys*9); S += bending_modes_surface(&xyz, N_mode, b + iCoordSys*N_mode, BM + iCoordSys*BM_N_SAMPLE*BM_N_SAMPLE, BM_N_SAMPLE, BM_radius); @ \end{itemize} Both above operations are combined together with: \begin{itemize} \item Zernike modes: <>= <> <> @ \item bending modes: <>= <> <> @ \end{itemize} The surface derivatives follow the same procedure than the surface but before adding the aspheric derivatives to the conic derivatives, the conic gradient vector is tranformed back into the motion CS: \begin{itemize} \item Zernike modes: <>= // CONIC >>> backward_transform_centered(&xyz, &xyz, d__conic_R+iCoordSys*9); K = partial_x_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x)); L = partial_y_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x)); M = partial_z_conic_surface(); // CONIC <<< forward_transform_centered(&xyz, &xyz, d__conic_R+iCoordSys*9); Snormal.x = K; Snormal.y = L; Snormal.z = M; forward_transform_centered(&Snormal, &Snormal, d__conic_R+iCoordSys*9); K = Snormal.x; L = Snormal.y; M = Snormal.z; K += partial_x_zernike_surface(&xyz, &d__origin, max_n, d__cx + iCoordSys*b_n_mode, R); L += partial_y_zernike_surface(&xyz, &d__origin, max_n, d__cy + iCoordSys*b_n_mode, R); M += partial_z_zernike_surface(); @ \item bending modes: <>= // CONIC >>> backward_transform_centered(&xyz, &xyz, d__conic_R+iCoordSys*9); K = partial_x_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x)); L = partial_y_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x)); M = partial_z_conic_surface(); // CONIC <<< forward_transform_centered(&xyz, &xyz, d__conic_R+iCoordSys*9); Snormal.x = K; Snormal.y = L; Snormal.z = M; forward_transform_centered(&Snormal, &Snormal, d__conic_R+iCoordSys*9); K = Snormal.x; L = Snormal.y; M = Snormal.z; K += partial_x_bending_modes_surface(&xyz, N_mode, b + iCoordSys*N_mode, BM + iCoordSys*BM_N_SAMPLE*BM_N_SAMPLE, BM_N_SAMPLE, BM_radius); L += partial_y_bending_modes_surface(&xyz, N_mode, b + iCoordSys*N_mode, BM + iCoordSys*BM_N_SAMPLE*BM_N_SAMPLE, BM_N_SAMPLE, BM_radius); M += partial_z_bending_modes_surface(); @ \end{itemize} The surface and the surface derivatives are computed together with: \begin{itemize} \item Zernike modes: <>= <> K = partial_x_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x)); L = partial_y_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x)); M = partial_z_conic_surface(); <> Snormal.x = K; Snormal.y = L; Snormal.z = M; forward_transform_centered(&Snormal, &Snormal, d__conic_R+iCoordSys*9); K = Snormal.x; L = Snormal.y; M = Snormal.z; K += partial_x_zernike_surface(&xyz, &d__origin, max_n, d__cx + iCoordSys*b_n_mode, R); L += partial_y_zernike_surface(&xyz, &d__origin, max_n, d__cy + iCoordSys*b_n_mode, R); M += partial_z_zernike_surface(); @ \item bending modes: <>= <> K = partial_x_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x)); L = partial_y_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x)); M = partial_z_conic_surface(); <> Snormal.x = K; Snormal.y = L; Snormal.z = M; forward_transform_centered(&Snormal, &Snormal, d__conic_R+iCoordSys*9); K = Snormal.x; L = Snormal.y; M = Snormal.z; K += partial_x_bending_modes_surface(&xyz, N_mode, b + iCoordSys*N_mode, BM + iCoordSys*BM_N_SAMPLE*BM_N_SAMPLE, BM_N_SAMPLE, BM_radius); L += partial_y_bending_modes_surface(&xyz, N_mode, b + iCoordSys*N_mode, BM + iCoordSys*BM_N_SAMPLE*BM_N_SAMPLE, BM_N_SAMPLE, BM_radius); M += partial_z_bending_modes_surface(); @ \end{itemize} This is the algorithm to compute the coordinates of the intersection of a ray with a conic surface: <>= k = klm.x; l = klm.y; m = klm.z; if (m==0) { return; } s0 = -xyz.z/m; x1 = xyz.x + k*s0; y1 = xyz.y + l*s0; d__ray[ij].optical_path_length = s0; s0 = s1 = 0; for (j=0; j> <> dSds = K*k + L*l + M*m; if (dSds==0) { break; } s1 = s0 - S/dSds; if (abs(s1-s0)> d__ray[ij].optical_path_length += s1; d__ray[ij].n_iteration = j; break; } s0 = s1; } @ <>= S = conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x)); @ <>= K = partial_x_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x)); L = partial_y_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x)); M = partial_z_conic_surface(); @ <>= k = klm_GS.x; l = klm_GS.y; m = klm_GS.z; //if (m==0) { return; } //s0 = -xyz_GS.z/m; x1 = xyz_GS.x;// + k*s0; y1 = xyz_GS.y;// + l*s0; z1 = xyz_GS.z; //d__ray[ij].optical_path_length += s0; s0 = s1 = 0; for (j=0; j> <> dSds = K*k + L*l + M*m; if (dSds==0) s1 = s0; else s1 = s0 - S/dSds; if (abs(s1-s0)> K += Snormal.x; L += Snormal.y; M += Snormal.z; d__ray[ij].optical_path_length += s1; d__ray[ij].n_iteration = j; break; } s0 = s1; } @ <>= S = zernike_surface(&xyz_GS, &d__origin, max_n, d__a + iCoordSys*n_mode, R); @ <>= K = partial_x_zernike_surface(&xyz_GS, &d__origin, max_n, d__cx + iCoordSys*b_n_mode, R); L = partial_y_zernike_surface(&xyz_GS, &d__origin, max_n, d__cy + iCoordSys*b_n_mode, R); M = partial_z_zernike_surface(); @ \subsection{Ray tracing without vignetting} \label{sec:ray-tracing-no-vignet} \index{gmtMirrors!gmt\_m1!traceall} <>= void gmt_m1::traceall(bundle *rays) { rays->V.area = area0*rays->N_BUNDLE*area_fraction; dim3 blockDim, gridDim; rtd R2, Rri2; R2 = D_clear*D_clear*0.25; Rri2 = R2*ri*ri*0; <> } @ where <>= /* if (V!=NULL) { */ /* blockDim = dim3(N_THREAD,N_THREAD); */ /* gridDim = dim3(rays->N_RAY/N_THREAD2+1,1,rays->N_BUNDLE); */ /* intersection LLL gridDim , blockDim RRR (V->m, rays->d__ray, rays->N_RAY); */ /* } */ V = &(rays->V); // fprintf(stdout,"R2=%5.2f - Rri2=%5.2f\n",R2,Rri2); blockDim = dim3(1,1); gridDim = dim3(1,1,rays->N_BUNDLE); m1_trace_chief_kernel LLL gridDim , blockDim RRR (V->m, rays->d__chief_ray, 1, conic_CS.d__R, conic_CS.d__origin, conic_k, conic_c, d__conic_origin, rigid_body_CS.d__R, rigid_body_CS.d__origin, motion_CS.d__R, motion_CS.d__origin); blockDim = dim3(N_THREAD,N_THREAD); gridDim = dim3(rays->N_RAY/N_THREAD2+1,1,rays->N_BUNDLE); /*HANDLE_ERROR( cudaMemset(V->m, 0, sizeof(char)*rays->N_RAY*rays->N_BUNDLE ) ); HANDLE_ERROR( cudaMemset(rays->d__piston_mask, 0, sizeof(int)*rays->N_RAY*rays->N_BUNDLE ) ); */ if (ZS!=NULL) m1_trace_all_kernel LLL gridDim , blockDim RRR (V->m, rays->d__ray, rays->N_RAY, Rri2, R2, aperture_CS.d__R, aperture_CS.d__origin, conic_CS.d__R, conic_CS.d__origin, d__conic_k, d__conic_c, d__conic_origin, rigid_body_CS.d__R, rigid_body_CS.d__origin, motion_CS.d__R, motion_CS.d__origin, ZS->max_n, ZS->n_mode, ZS->d__a, ZS->d__cx, ZS->d__cy, 0.5*D_full, rays->d__chief_ray, rays->d__piston_mask, d__valid_segments, d__segment_reflectivity); /* if (BS!=NULL) m1_bm_trace_kernel LLL gridDim , blockDim RRR (V->m, rays->d__ray, rays->N_RAY, Rri2, R2, aperture_CS.d__R, aperture_CS.d__origin, conic_CS.d__R, conic_CS.d__origin, d__conic_k, d__conic_c, d__conic_origin, rigid_body_CS.d__R, rigid_body_CS.d__origin, motion_CS.d__R, motion_CS.d__origin, BS.n_mode, BS.d__b, BS.d__BMS, BS.BM_N_SAMPLE, BS.BM_radius, rays->d__chief_ray, rays->d__piston_mask, d__valid_segments, d__segment_reflectivity); intersection LLL gridDim , blockDim RRR (V->m, rays->d__ray, rays->N_RAY); float previous_nnz = rays->V.nnz; //fprintf(stdout,"M%d previous nnz rays: %f\n",M_ID,previous_nnz); rays->V.set_filter_quiet(); //fprintf(stdout,"M%d current nnz rays: %f\n",M_ID,rays->V.nnz); if (previous_nnzN_RAY_TOTAL) rays->V.area *= rays->V.nnz/previous_nnz; //fprintf(stdout,"M%d current area: %f\n",M_ID,rays->V.area); */ @ and the kernel <>= __global__ void m1_trace_all_kernel(char *mask, ray *d__ray, int N_RAY, rtd inner2, rtd outer2, rtd *d__aperture_R, vector *d__aperture_origin, rtd *d__conic_R, vector *d__conic_origin, rtd *d__Fk, rtd *d__Fc, vector *d__conic_self_origin, rtd *d__rigid_body_R, vector *d__rigid_body_origin, rtd *d__motion_R, vector *d__motion_origin, int max_n, int n_mode, rtd *d__a, rtd *d__cx, rtd *d__cy, rtd R, ray *d__chief_ray, int *d__piston_mask, char *d__valid_segments, float *reflectivity) { int i, j, ij, iCoordSys, iSource, b_n_mode; rtd rho2; rtd x, y, z, x1, y1, s0, k, l, m; rtd s1, S, K, L ,M, dSds; rtd G2, a; rtd Fc, Fk; vector d__origin, xyz, klm, xyz_GS, klm_GS, Snormal; i = blockIdx.x * blockDim.x + threadIdx.x; j = threadIdx.y; // iCoordSys = blockIdx.y; ij = j * gridDim.x * blockDim.x + i; iSource = blockIdx.z; j = ij; ij += iSource*N_RAY; b_n_mode = max_n*(max_n+1)*0.5; if ( ( j0) { d__origin = d__conic_self_origin[iCoordSys]; <> { mask[ij] = 1; d__piston_mask[ij] = iCoordSys + 1; <> d__ray[ij].optical_path_difference += d__ray[ij].optical_path_length - d__chief_ray[iSource].optical_path_length; d__ray[ij].throughput *= reflectivity[iCoordSys]; return; } } } } } @ \subsection{Segment selection} \label{sec:segment-selection} Segments can be removed from the M1 or M2 models with the routine: \index{gmtMirrors!gmt\_m1!remove} <>= void gmt_m1::remove(int *seg_ID, int N_ID) { <> } @ where [[seg_ID]] is an array of length [[N_ID]] with the ID number of the segments to deactivate, and <>= int k; char valid[7]; area_fraction = (7-N_ID)/7.0; memset(valid, 1, sizeof(char)*7); for (k=0;k>= void gmt_m1::keep(int *seg_ID, int N_ID) { <> } @ where [[seg_ID]] is an array of length [[N_ID]] with the ID number of the segments to keep, and <>= int k; char valid[7]; area_fraction = N_ID/7.0; memset(valid, 0, sizeof(char)*7); for (k=0;k>= void gmt_m1::track(float *d__x, float *d__y, float *d__z, int N, int idx) { <> } @ with <>= dim3 blockDim(N_THREAD,1); dim3 gridDim(N/N_THREAD+1,1); track_kernel LLL gridDim, blockDim RRR (d__x, d__y, d__z, N, idx, rigid_body_CS.d__R, rigid_body_CS.d__origin, motion_CS.d__R, motion_CS.d__origin); @ and <>= __global__ void track_kernel(float *x, float *y, float *z, const int N, const int iCoordSys, rtd *d__rigid_body_R, vector *d__rigid_body_origin, rtd *d__motion_R, vector *d__motion_origin) { int i; vector v; i = blockIdx.x * blockDim.x + threadIdx.x; if (i>= void gmt_m1::locate(float *d__x, float *d__y, float *d__z, int N, int idx) { <> } @ with <>= dim3 blockDim(N_THREAD,1); dim3 gridDim(N/N_THREAD+1,1); locate_kernel LLL gridDim, blockDim RRR (d__x, d__y, d__z, N, idx, rigid_body_CS.d__R, rigid_body_CS.d__origin, motion_CS.d__R, motion_CS.d__origin); @ and <>= __global__ void locate_kernel(float *x, float *y, float *z, const int N, const int iCoordSys, rtd *d__rigid_body_R, vector *d__rigid_body_origin, rtd *d__motion_R, vector *d__motion_origin) { int i; vector v; i = blockIdx.x * blockDim.x + threadIdx.x; if (i>> forward_transform(&v, &v, d__rigid_body_R+iCoordSys*9, d__rigid_body_origin+iCoordSys); // MOTION >>> forward_transform(&v, &v, d__motion_R+iCoordSys*9, d__motion_origin+iCoordSys); x[i] = v.x; y[i] = v.y; z[i] = v.z; } } @ \subsection{Global tip--tilt} \label{sec:global-tip-tilt} \def\ijk{\ensuremath{\left[\vec i,\vec j, \vec k\right]}} \def\uvw{\ensuremath{\left[\vec u,\vec v, \vec w\right]_k}} \def\uvwp{\ensuremath{\left[\vec u^\prime,\vec v^\prime, \vec w^\prime\right]}} \def\RS{\ensuremath{R_{{\cal S},k}}} \def\RT{\ensuremath{R_{\cal T}}} \def\RM{\ensuremath{R_{\cal M}}} The GCS is defined with the 3 ortho--normal unit vectors \ijk. The segment coordinate systems also called rigid body CSs are defined with the 3 ortho-normal vectors \uvw. \uvw are derived from \ijk with the matrix \RS \begin{equation} \label{eq:3} \uvw = \RS^T \ijk. \end{equation} Within the rigid body CSs, motion CSs (\uvwp) are defined with \begin{equation} \label{eq:4} \uvwp = \RM^T \uvw = \RM^T \RS^T \ijk, \end{equation} where \RM is the 3D rotation matrix of each segment. When a global tip--tilt is applied to M1, \ijk becomes $\RT^T\ijk$ with \RT the tip--tilt matrix, and the motion CSs are given by \begin{equation} \label{eq:5} \uvwp = \RS^T\RT^T\ijk. \end{equation} From Eq.~(\ref{eq:4}) and Eq.~(\ref{eq:5}), the 3D rotation matrix of each segment as a result of a global tip--tilt is given by \begin{equation} \label{eq:6} \RM = \RS^T \RT \RS. \end{equation} The transpose of Eq.~(\ref{eq:6}) is implemented below, remembering that the transpose of the rotation matrices are stored into CEO. <>= cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, 3,3,3, &alpha, TT_CS.d__R, 3, rigid_body_CS.d__R+k9, 3, &beta, d__C, 3); cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, 3,3,3, &alpha, rigid_body_CS.d__R+k9, 3, d__C, 3, &beta, d__D, 3); cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, 3,3,3, &alpha, motion_CS.d__R+k9, 3, d__D, 3, &beta, d__C, 3); HANDLE_ERROR( cudaMemcpy( motion_CS.d__R+k9, d__C, sizeof(double)*9, cudaMemcpyDeviceToDevice ) ); HANDLE_ERROR( cudaMemcpy( motion_CS.R+k9, motion_CS.d__R+k9, sizeof(double)*9, cudaMemcpyDeviceToHost ) ); /* int i, j, idx = 0; for (i=0;i<3;i++) { fprintf(stdout,"|| "); for (j=0;j<3;j++) { idx = j + i*3 + k9; fprintf(stdout,"%+.4e ",motion_CS.R[idx]); } fprintf(stdout,"||\n"); } */ @ From \RM, the Euler angles corresponding to the 3 rotations of each segment are derived \begin{eqnarray} \label{eq:7} \alpha &=& \arctan\left( \RM(32) \over \RM(33) \right) \\\nonumber \beta &=& \arcsin\left( -\RM(31) \right) \\\nonumber \gamma &=& \arctan\left( \RM(21) \over \RM(11) \right) \end{eqnarray} <>= motion_CS.euler_angles[k].x = atan2(motion_CS.R[7+k9],motion_CS.R[8+k9]); motion_CS.euler_angles[k].y = asin(-1.0*motion_CS.R[6+k9]); motion_CS.euler_angles[k].z = atan2(motion_CS.R[3+k9],motion_CS.R[k9]); @ The tip--tilt transformation is also applied to the origin of each segment $\vec O_k$ and the difference between the two origins is the translation of the motion CS, i.e. \begin{equation} \label{eq:8} \vec O_k^\prime = \vec O_k - \vec O_R - \RT^T\left( \vec O_k -\vec O_R \right) \end{equation} where $\vec O_R$ is the rotation point. <>= vector o; forward_transform(&o, rigid_body_CS.origin + k, TT_CS.R, TT_CS.origin); o.x = rigid_body_CS.origin[k].x - TT_CS.origin->x - o.x; o.y = rigid_body_CS.origin[k].y - TT_CS.origin->y - o.y; o.z = rigid_body_CS.origin[k].z - TT_CS.origin->z - o.z; @ The new segment origin is transformed into the motion CS and added to the current motion CS locations: <>= forward_transform_centered(&o, &o, rigid_body_CS.R + k9); motion_CS.origin[k].x += o.x; motion_CS.origin[k].y += o.y; motion_CS.origin[k].z += o.z; HANDLE_ERROR( cudaMemcpy( motion_CS.d__origin+k, motion_CS.origin+k, sizeof(vector), cudaMemcpyHostToDevice ) ); /* fprintf(stdout,"#%d >> MOTION CS origins [micron] : %+.3e, %+.3e, %+.3e\n",k, motion_CS.origin[k].x*1e6, motion_CS.origin[k].y*1e6, motion_CS.origin[k].z*1e6); float r2d; r2d = 1000*3600*180.0/PI; fprintf(stdout,"#%d >> MOTION CS Euler angles [mas]: %+.3e, %+.3e, %+.3e\n",k, r2d*motion_CS.euler_angles[k].x, r2d*motion_CS.euler_angles[k].y, r2d*motion_CS.euler_angles[k].z); */ @ \index{gmtMirrors!gmt\_m1!global\_tiptilt} <>= void gmt_m1::global_tiptilt(float tip, float tilt) { <> } @ with <>= sprintf(tag,"GLOBAL TIP-TILT"); vector TT_CS_origin, TT_CS_euler_angles; TT_CS_origin.x = TT_CS_origin.y = TT_CS_origin.z = 0.0; TT_CS_euler_angles.x = TT_CS_euler_angles.y = TT_CS_euler_angles.z = 0.0; TT_CS.setup(&TT_CS_origin, &TT_CS_euler_angles, 1, tag); cublasCreate(&handle); HANDLE_ERROR( cudaMalloc((void**)&d__C, sizeof(double)*9 ) ); HANDLE_ERROR( cudaMalloc((void**)&d__D, sizeof(double)*9 ) ); @ <>= TT_CS.cleanup(); cublasDestroy(handle); HANDLE_ERROR( cudaFree( d__C ) ); HANDLE_ERROR( cudaFree( d__D ) ); @ and <>= double alpha, beta; alpha = 1.0; beta = 0.0; int k, k9; vector origin, euler_angles, o; origin.x = origin.y = origin.z = 0.0; euler_angles.x = tip; euler_angles.y = tilt; euler_angles.z = 0.0; TT_CS.update(origin,euler_angles,0); for (k=0;k> } @ \subsection{Segment reflectivity} \label{sec:seg-reflect} The segment reflectivity is set with: \index{gmtMirrors!gmt\_m1!set\_reflectivity} <>= void gmt_m1::set_reflectivity(float *reflectivity) { <> } @ with <>= int n_byte; n_byte = 7*sizeof(float); HANDLE_ERROR( cudaMemcpy( d__segment_reflectivity, reflectivity, n_byte, cudaMemcpyHostToDevice ) ); @ \section{GMT M2} \label{sec:gmt-m2} \index{gmtMirrors!gmt\_m2} A new structure to hold GMT M1 parameters and functions is defined. <>= struct gmt_m2 { <> <> void pointing_neutral(float tip, float tilt); void coma_neutral(float tip, float tilt); }; @ \subsection{Setup \& Cleanup} \label{sec:setup--cleanup-2} \index{gmtMirrors!gmt\_m2!setup} <>= void gmt_m2::setup(void) { // BS = NULL; ZS = NULL; <> <> } @ <>= void gmt_m2::setup(modes *_BS_) { BS = _BS_; ZS = NULL; <> <> } @ <>= void gmt_m2::setup(char *_filename_, int _N_, int _n_mode_) { BS.setup(_filename_, _N_, _n_mode_); ZS = NULL; <> <> } @ with <>= M_ID = 2; D_assembly = 3.168; D_full = 1.0425*1.005; // slightly larger to always encompass M1; needed for the Zernike D_clear = 1.0415; ri = 0.0; beta = -14.777498*PI/180.0; L = 1.08774; N = 7; conic_c = -1.0/4.1639009; conic_k = 1-0.71692784; height = 20.26247614; idx_offset = 3; @ \index{gmtMirrors!gmt\_m2!setup} <>= void gmt_m2::setup(zernikeS *_ZS_) { ZS = _ZS_; //BS = NULL; M_ID = 2; <> <> } @ Memory is freed with \index{gmtMirrors!gmt\_m2!cleanup} <>= void gmt_m2::cleanup(void) { INFO("@(CEO)>gmt_m2: freeing memory!\n"); <> } @ <>= void gmt_m2::update_conic_c(rtd *_conic_c_) { HANDLE_ERROR( cudaMemcpy( d__conic_c, _conic_c_, sizeof(rtd)*7, cudaMemcpyHostToDevice ) ); } @ <>= void gmt_m2::update_conic_k(rtd *_conic_k_) { HANDLE_ERROR( cudaMemcpy( d__conic_k, _conic_k_, sizeof(rtd)*7, cudaMemcpyHostToDevice ) ); } @ \subsection{Blocking} \label{sec:blocking-2} The rays blocking with M2 is computed with \index{gmtMirrors!gmt\_m2!blocking} <>= void gmt_m2::blocking(bundle *rays) { <> } @ \subsection{Ray tracing} \label{sec:ray-tracing-2} The rays propagation through M1 is computed with \index{gmtMirrors!gmt\_m2!trace} <>= void gmt_m2::trace(bundle *rays) { dim3 blockDim, gridDim; rtd R2, Rri2; R2 = D_clear*D_clear*0.25; Rri2 = R2*ri*ri; <> } void gmt_m2::traceall(bundle *rays) { rays->V.area = area0*rays->N_BUNDLE*area_fraction; dim3 blockDim, gridDim; rtd R2, Rri2; R2 = D_clear*D_clear*0.25; Rri2 = R2*ri*ri*0; <> } @ The rigid body motion parameters are updated with: \index{gmtMirrors!gmt\_m2!update} <>= void gmt_m2::update(vector _origin_, vector _euler_angles_,int idx) { <> } @ The rigid body motion are reset to 0 with: \index{gmtMirrors!gmt\_m2!reset} <>= void gmt_m2::reset(void) { <> } @ \subsection{Segment selection} \label{sec:segment-selection} Segments can be removed from the M1 or M2 models with the routine: \index{gmtMirrors!gmt\_m2!remove} <>= void gmt_m2::remove(int *seg_ID, int N_ID) { <> } @ One can also specifies the segments to keep \index{gmtMirrors!gmt\_m2!keep} <>= void gmt_m2::keep(int *seg_ID, int N_ID) { <> } @ \subsection{Mirror location tracker} \label{sec:mirr-locat-track} The [[track]] function transforms the coordinates $(x,y,z)$ in the motion CS of segment [[idx]] to the GCS. \index{gmtMirrors!gmt\_m2!track} <>= void gmt_m2::track(float *d__x, float *d__y, float *d__z, int N, int idx) { <> } @ The [[locate]] function transforms the coordinates $(x,y,z)$ in the GCS to the motion CS of segment [[idx]] \index{gmtMirrors!gmt\_m2:locate} <>= void gmt_m2::locate(float *d__x, float *d__y, float *d__z, int N, int idx) { <> } @ \subsection{Global tip--tilt} \label{sec:global-tip-tilt-1} \index{gmtMirrors!gmt\_m2!global\_tiptilt} <>= void gmt_m2::global_tiptilt(float tip, float tilt) { <> } @ \subsubsection{Tip--tilt neutral} \label{sec:tip-tilt-neutral} \index{gmtMirrors!gmt\_m2!pointing\_neutral} <>= void gmt_m2::pointing_neutral(float tip, float tilt) { double alpha, beta, *d__C; vector origin, euler_angles; cublasHandle_t handle; coordinate_system TT_CS; origin.x = origin.y = 0.0; origin.z = height - 4390.312E-03; euler_angles.x = tip; euler_angles.y = tilt; euler_angles.z = 0.0; TT_CS.setup(origin, euler_angles); cublasCreate(&handle); alpha = 1; beta = 0; HANDLE_ERROR( cudaMalloc((void**)&d__C, sizeof(double)*9 ) ); int k, k9; for (k=0;k> } cublasDestroy(handle); TT_CS.cleanup(); HANDLE_ERROR( cudaFree( d__C ) ); } @ \subsubsection{Coma neutral} \label{sec:coma-neutral} \index{gmtMirrors!gmt\_m2!coma\_neutral} <>= void gmt_m2::coma_neutral(float tip, float tilt) { double alpha, beta, *d__C; vector origin, euler_angles; cublasHandle_t handle; coordinate_system TT_CS; origin.x = origin.y = 0.0; origin.z = height - 2246.410E-03; euler_angles.x = tip; euler_angles.y = tilt; euler_angles.z = 0.0; TT_CS.setup(origin, euler_angles); cublasCreate(&handle); alpha = 1; beta = 0; HANDLE_ERROR( cudaMalloc((void**)&d__C, sizeof(double)*9 ) ); int k, k9; for (k=0;k> } cublasDestroy(handle); TT_CS.cleanup(); HANDLE_ERROR( cudaFree( d__C ) ); } @ \subsection{Segment reflectivity} \label{sec:seg-reflect} The segment reflectivity is set with: \index{gmtMirrors!gmt\_m2!set\_reflectivity} <>= void gmt_m2::set_reflectivity(float *reflectivity) { <> } @ \section{Bending modes} \label{sec:bending-modes} This section defines the routine to load the bending modes and interpolated the bending modes; The parameters and routines associated to the bending modes are gatherered in the \emph{bending\_modes} structure: \index{gmtMirrors!bending\_modes} <>= struct bending_modes { <> void setup(int n_mode); void setup(int n_mode, int N_SET); void setupKL(int n_mode); void setupKL(int n_mode, int N_SET); void setupPolish(int n_mode, int N_SET); void cleanup(void); void load(void); void load_reg(void); void load_KL(void); void load_polish(void); void nearest_neighbor(rtd *d__BMi, rtd *d__partial_x_BMi, rtd *d__partial_y_BMi, int NI, rtd di, int k_mode); void bilinear(rtd *d__BMi, rtd *d__partial_x_BMi, rtd *d__partial_y_BMi, int NI, rtd di, int k_mode); void update(rtd *b); }; @ The bending modes parameters are: \begin{itemize} \item the bending modes x and y coordinates <>= double *d__x_BM, *d__y_BM; @ \item the bending modes <>= double *d__BM; @ \item the bending modes surface <>= double *d__BMS; @ \item the radius of the circle circumscribing the bending modes <>= double BM_radius; @ \item the number of sampling points i.e. the length of [[x_BM]] and [[y_BM]] arrays <>= int BM_N_SAMPLE; @ \item the bending modes data ([[d__x_BM]],[[d__y_BM]] and [[d__BM]]) buffer <>= double *d__BM_buffer; @ \item the number of bending modes <>= int n_mode; @ \item the bending modes coefficients [[b]] <>= rtd *b, *d__b; @ \item the number of bending modes set [[N]] <>= int N; @ \end{itemize} \subsection{Setup \& Cleanup} \label{sec:setup--cleanup} The bending modes are loaded with: \index{gmtMirrors!bending\_modes!setup} <>= void bending_modes::setup(int _n_mode_) { N = 1; n_mode = _n_mode_; <> load_reg(); <> } @ or for a set of [[N]] bending modes <>= void bending_modes::setup(int _n_mode_, int N_SET) { N = N_SET; n_mode = _n_mode_; <> load_reg(); <> } @ \index{gmtMirrors!bending\_modes!setupKL} <>= void bending_modes::setupKL(int _n_mode_) { N = 1; n_mode = _n_mode_; <> load_KL(); <> } @ or for a set of [[N]] bending modes <>= void bending_modes::setupKL(int _n_mode_, int N_SET) { N = N_SET; n_mode = _n_mode_; <> load_KL(); <> } @ or for the polishing map <>= void bending_modes::setupPolish(int _n_mode_, int N_SET) { N = N_SET; n_mode = _n_mode_; <> load_polish(); <> } @ with <>= int n_byte = sizeof(rtd)*n_mode*N; b = (rtd *)malloc( n_byte ); HANDLE_ERROR( cudaMalloc((void**)&d__b, n_byte ) ); HANDLE_ERROR( cudaMemcpy( d__b, b, n_byte, cudaMemcpyHostToDevice ) ); HANDLE_ERROR( cudaMemset(d__b, 0, n_byte ) ); <>= n_byte = BM_N_SAMPLE*BM_N_SAMPLE*N*sizeof(double); HANDLE_ERROR( cudaMalloc((void**)&d__BMS, n_byte ) ); HANDLE_ERROR( cudaMemset(d__BMS, 0, n_byte ) ); @ Memory is freed with \index{gmtMirrors!bending\_modes!cleanup} <>= void bending_modes::<> <>= cleanup(void) { INFO("@(CEO)>bending_modes: freeing memory!\n"); free( b ); HANDLE_ERROR( cudaFree( d__b ) ); HANDLE_ERROR( cudaFree( d__BM_buffer ) ); HANDLE_ERROR( cudaFree( d__BMS ) ); } @ \subsection{Update} \label{sec:update} The bending mode coefficients [[d__b]] and the associated surface [[d__BMS]] are updated with: <>= void bending_modes::update(rtd *b) { int n_byte = sizeof(rtd)*n_mode*N; HANDLE_ERROR( cudaMemcpy( d__b, b, n_byte, cudaMemcpyHostToDevice ) ); dim3 blockDim(N_THREAD,N_THREAD); dim3 gridDim(BM_N_SAMPLE/N_THREAD+1,BM_N_SAMPLE/N_THREAD+1,N); surface_update LLL gridDim, blockDim RRR (d__BMS, n_mode, d__b, d__BM, BM_N_SAMPLE); } @ with the kernel <>= __global__ void surface_update(double *BMS, const int N_mode, rtd *b, double *BM, int N_SAMPLE) { int i, j, k, kz, l, k_mode, o; rtd *z; i = blockIdx.x * blockDim.x + threadIdx.x; j = blockIdx.y * blockDim.y + threadIdx.y; l = blockIdx.z; if ( (i>= void bending_modes::load(void) { FILE * pFile; long lSize; double * buffer; size_t result; char path[256]; sprintf(path, "%sgmtMirrors/bendingModes.bin",getenv("CEOPATH")); INFO("@(CEO)>bending_modes: loading bending modes from %s\n",path); pFile = fopen ( path , "rb" ); if (pFile==NULL) {fputs ("File error",stderr); exit (1);} // obtain file size: fseek (pFile , 0 , SEEK_END); lSize = ftell (pFile)/sizeof(double); rewind (pFile); // fprintf(stdout,"lSize=%ld\n",lSize); // allocate memory to contain the whole file: buffer = (double*) malloc (lSize*sizeof(double)); if (buffer == NULL) {fputs ("Memory error",stderr); exit (2);} // copy the file into the buffer: result = fread (buffer,sizeof(double),lSize,pFile); // fprintf(stdout,"result=%ld\n",result); if (result != lSize) {fputs ("Reading error",stderr); exit (3);} /* the whole file is now loaded in the memory buffer. */ BM_radius = 4.181; BM_N_SAMPLE = 27685; HANDLE_ERROR( cudaMalloc((void**)&d__BM_buffer, lSize*sizeof(double)) ); HANDLE_ERROR( cudaMemcpy( d__BM_buffer , buffer, lSize*sizeof(double), cudaMemcpyHostToDevice ) ); d__x_BM = d__BM_buffer; d__y_BM = d__BM_buffer + BM_N_SAMPLE; d__BM = d__BM_buffer + BM_N_SAMPLE*2; // terminate fclose (pFile); free (buffer); } @ \subsubsection{Regular grid} \label{sec:regular-grid} The bending modes are interpolated on a regular grid are saved in the file \emph{bendingModesReg.bin}. The regular grid is made of $169\times 169$ points. The range of the grid is $\pm$4.231m. The bending modes are loaded into the \emph{gmt\_m1} structure with: \index{gmtMirrors!bending\_modes!load\_reg} <>= void bending_modes::load_reg(void) { FILE * pFile; long lSize; double * buffer; size_t result; char path[256]; sprintf(path, "%sgmtMirrors/bendingModesReg.bin",getenv("CEOPATH")); INFO("@(CEO)>bending_modes: loading bending modes from %s\n",path); pFile = fopen ( path , "rb" ); if (pFile==NULL) {fputs ("File error",stderr); exit (1);} // obtain file size: fseek (pFile , 0 , SEEK_END); lSize = ftell (pFile)/sizeof(double); rewind (pFile); // fprintf(stdout,"lSize=%ld\n",lSize); // allocate memory to contain the whole file: buffer = (double*) malloc (lSize*sizeof(double)); if (buffer == NULL) {fputs ("Memory error",stderr); exit (2);} // copy the file into the buffer: result = fread (buffer,sizeof(double),lSize,pFile); // fprintf(stdout,"result=%ld\n",result); if (result != lSize) {fputs ("Reading error",stderr); exit (3);} /* the whole file is now loaded in the memory buffer. */ BM_radius = 0.5*8.462;//4.181; BM_N_SAMPLE = 169; HANDLE_ERROR( cudaMalloc((void**)&d__BM_buffer, lSize*sizeof(double)) ); HANDLE_ERROR( cudaMemcpy( d__BM_buffer , buffer, lSize*sizeof(double), cudaMemcpyHostToDevice ) ); d__BM = d__BM_buffer; // terminate fclose (pFile); free (buffer); } @ \subsubsection{Polishing map} \label{sec:polishing-map} The polishing map is interpolated on a regular grid are saved in the file \emph{polishingMap.bin}. The regular grid is made of $483\times 483$ points. The range of the grid is $\pm$4.26m. The bending modes are loaded into the \emph{gmt\_m1} structure with: \index{gmtMirrors!bending\_modes!load\_polish} <>= void bending_modes::load_polish(void) { FILE * pFile; long lSize; double * buffer; size_t result; char path[256]; sprintf(path, "%sgmtMirrors/polishingMap.bin",getenv("CEOPATH")); INFO("@(CEO)>bending_modes: loading bending modes from %s\n",path); pFile = fopen ( path , "rb" ); if (pFile==NULL) {fputs ("File error",stderr); exit (1);} // obtain file size: fseek (pFile , 0 , SEEK_END); lSize = ftell (pFile)/sizeof(double); rewind (pFile); // fprintf(stdout,"lSize=%ld\n",lSize); // allocate memory to contain the whole file: buffer = (double*) malloc (lSize*sizeof(double)); if (buffer == NULL) {fputs ("Memory error",stderr); exit (2);} // copy the file into the buffer: result = fread (buffer,sizeof(double),lSize,pFile); // fprintf(stdout,"result=%ld\n",result); if (result != lSize) {fputs ("Reading error",stderr); exit (3);} /* the whole file is now loaded in the memory buffer. */ BM_radius = 0.5*8.52; BM_N_SAMPLE = 483; HANDLE_ERROR( cudaMalloc((void**)&d__BM_buffer, lSize*sizeof(double)) ); HANDLE_ERROR( cudaMemcpy( d__BM_buffer , buffer, lSize*sizeof(double), cudaMemcpyHostToDevice ) ); d__BM = d__BM_buffer; // terminate fclose (pFile); free (buffer); } @ \subsubsection{Karhunen--Loeve} \label{sec:regular-grid} The Karhunen--Loeve modes are interpolated on a regular grid are saved in the file \emph{KarhunenLoeveModes.bin}. The regular grid is made of $256\times 256$ points. The range of the grid is $\pm$4.25m. The bending modes are loaded into the \emph{gmt\_m2} structure with: \index{gmtMirrors!bending\_modes!load\_bending\_modes\_reg} <>= void bending_modes::load_KL(void) { FILE * pFile; long lSize; double * buffer; size_t result; char path[256]; sprintf(path, "%sgmtMirrors/KarhunenLoeveModes.bin",getenv("CEOPATH")); INFO("@(CEO)>bending_modes: loading bending modes from %s\n",path); pFile = fopen ( path , "rb" ); if (pFile==NULL) {fputs ("File error",stderr); exit (1);} // obtain file size: fseek (pFile , 0 , SEEK_END); lSize = ftell (pFile)/sizeof(double); rewind (pFile); // fprintf(stdout,"lSize=%ld\n",lSize); // allocate memory to contain the whole file: buffer = (double*) malloc (lSize*sizeof(double)); if (buffer == NULL) {fputs ("Memory error",stderr); exit (2);} // copy the file into the buffer: result = fread (buffer,sizeof(double),lSize,pFile); // fprintf(stdout,"result=%ld\n",result); if (result != lSize) {fputs ("Reading error",stderr); exit (3);} /* the whole file is now loaded in the memory buffer. */ BM_radius = 0.5*1.0425*1.1;//4.25; BM_N_SAMPLE = 256; HANDLE_ERROR( cudaMalloc((void**)&d__BM_buffer, lSize*sizeof(double)) ); HANDLE_ERROR( cudaMemcpy( d__BM_buffer , buffer, lSize*sizeof(double), cudaMemcpyHostToDevice ) ); d__BM = d__BM_buffer; // terminate fclose (pFile); free (buffer); } @ \subsection{Nearest neighbor interpolation} \label{sec:near-neighb-interp} The bending modes are interpolation on a regular $[[NI]]\times[[NI]]$ grid with spacing [[di]]. The bending modes are computed on an irregular grid within a circle of diameter 8.362m. If a point lies outside the circle, the interpolation returns 0. The routine also returns the x and y derivatives of the bending modes using the centered finite difference: \begin{equation} \label{eq:17} {\partial B(x,y) \over \partial x} = {B(x+h,y) - B(x-h,y) \over 2h}, \end{equation} with $h=5cm$ (the approximate spacing of the irregular grid). The interpolation algorithm uses a simple nearest neighbor interpolation method where the interpolated value at a given point is given by the value of the bending modes of the nearest point on the origin mesh. \index{gmtMirrors!bending\_modes!nearest\_neighbor} <>= void bending_modes::nearest_neighbor(rtd *d__BMi, rtd *d__partial_x_BMi, rtd *d__partial_y_BMi, int NI, rtd di, int k_mode) { dim3 blockDim(N_THREAD,N_THREAD); dim3 gridDim(NI/N_THREAD+1,NI/N_THREAD+1); nearest_neighbor_kernel LLL gridDim, blockDim RRR (d__BMi, d__partial_x_BMi, d__partial_y_BMi, NI, di, d__x_BM, d__y_BM, d__BM, BM_N_SAMPLE, BM_radius, k_mode); } @ with the kernel <>= __global__ void nearest_neighbor_kernel(rtd *BMi, rtd *partial_x_BMi, rtd *partial_y_BMi, const int NI, const rtd di, double *x_BM, double *y_BM, double *BM, const int BM_N_SAMPLE, const double BM_radius, const int k_mode) { int i, j, k, k_BM, t, o, t_E, t_N, t_W, t_S; rtd xi, yi, rho_t, rho_k, rhoi, h, xi_E, yi_N, xi_W, yi_S, rho_k_E, rho_k_N, rho_k_W, rho_k_S, rho_t_E, rho_t_N, rho_t_W, rho_t_S; i = blockIdx.x * blockDim.x + threadIdx.x; j = blockIdx.y * blockDim.y + threadIdx.y; if ( (i> BMi[k] = 0; partial_x_BMi[k] = 0; partial_y_BMi[k] = 0; if (rhoi<=BM_radius) { <> o = (k_mode-1)*BM_N_SAMPLE; BMi[k] = BM[t+o]; partial_x_BMi[k] = BM[t_W+o] - BM[t_E+o]; partial_x_BMi[k] /= (x_BM[t_W] - x_BM[t_E]); partial_y_BMi[k] = BM[t_N+o] - BM[t_S+o]; partial_y_BMi[k] /= (y_BM[t_N] - y_BM[t_S]); } } } @ with the grid coordinates definition <>= xi = i - 0.5*(NI-1); yi = j - 0.5*(NI-1); xi *= di; yi *= di; rhoi = hypot(xi,yi); h = 4*5e-2; xi_E = xi + h; yi_N = yi + h; xi_W = xi - h; yi_S = yi - h; @ and with the nearest neighbor algorithm <>= t = 0; rho_t = hypot(x_BM[t]-xi , y_BM[t]-yi); t_E = 0; rho_t_E = hypot(x_BM[t_E]-xi_E, y_BM[t_E]-yi); t_N = 0; rho_t_N = hypot(x_BM[t_N]-xi , y_BM[t_N]-yi_N); t_W = 0; rho_t_W = hypot(x_BM[t_W]-xi_W, y_BM[t_W]-yi); t_S = 0; rho_t_S = hypot(x_BM[t_S]-xi , y_BM[t_S]-yi_S); for (k_BM=0; k_BM>= void bending_modes::bilinear(rtd *d__BMi, rtd *d__partial_x_BMi, rtd *d__partial_y_BMi, int NI, rtd di, int k_mode) { dim3 blockDim(N_THREAD,N_THREAD); dim3 gridDim(NI/N_THREAD+1,NI/N_THREAD+1); bilinear_kernel LLL gridDim, blockDim RRR (d__BMi, d__partial_x_BMi, d__partial_y_BMi, NI, di, d__BM, BM_N_SAMPLE, BM_radius, k_mode); } @ with the kernels <>= __global__ void bilinear_kernel(rtd *BMi, rtd *partial_x_BMi, rtd *partial_y_BMi, const int NI, const rtd di, double *BM, const int BM_N_SAMPLE, const double BM_radius, const int k_mode) { int i, j, k, ndx, o; rtd xi, yi, rhoi, s, t, scale, fs, ft, onemt, onems, h; rtd *z; i = blockIdx.x * blockDim.x + threadIdx.x; j = blockIdx.y * blockDim.y + threadIdx.y; if ( (i> <> <> } } } @ where <>= s = (scale*(i + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1); t = (scale*(j + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1); <> BMi[k] = ( z[ndx]*onemt + z[ndx+1]*t )*onems + ( z[ndx+BM_N_SAMPLE]*onemt + z[ndx+BM_N_SAMPLE+1]*t )*s; @ <>= s = (scale*(i + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1) + 1; t = (scale*(j + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1); <> partial_x_BMi[k] = ( z[ndx]*onemt + z[ndx+1]*t )*onems + ( z[ndx+BM_N_SAMPLE]*onemt + z[ndx+BM_N_SAMPLE+1]*t )*s; s = (scale*(i + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1) - 1; t = (scale*(j + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1); <> partial_x_BMi[k] -= ( z[ndx]*onemt + z[ndx+1]*t )*onems + ( z[ndx+BM_N_SAMPLE]*onemt + z[ndx+BM_N_SAMPLE+1]*t )*s; partial_x_BMi[k] /= 2*h; @ <>= s = (scale*(i + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1); t = (scale*(j + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1) + 1; <> partial_y_BMi[k] = ( z[ndx]*onemt + z[ndx+1]*t )*onems + ( z[ndx+BM_N_SAMPLE]*onemt + z[ndx+BM_N_SAMPLE+1]*t )*s; s = (scale*(i + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1); t = (scale*(j + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1) - 1; <> partial_y_BMi[k] -= ( z[ndx]*onemt + z[ndx+1]*t )*onems + ( z[ndx+BM_N_SAMPLE]*onemt + z[ndx+BM_N_SAMPLE+1]*t )*s; partial_y_BMi[k] /= 2*h; @ and with <>= if (s<0) { s=0.0; } fs = floorf(s); if (t<0) { t=0.0; } ft = floorf(t); ndx = __float2int_rd( ft + fs*BM_N_SAMPLE ); if (s==(BM_N_SAMPLE-1)) { fs -= 1; ndx -= BM_N_SAMPLE; } if (t==(BM_N_SAMPLE-1)) { ft -= 1; ndx -= 1; } s -= fs; t -= ft; onemt = 1 - t; onems = 1 - s; @ \subsubsection{Bending mode surface} \label{sec:bending-mode-surface} A bending mode surface $S(\vec r)$ is given by a linear combination of bending modes $B_i(\vec r)$: \begin{equation} \label{eq:21} S(\vec r) = \sum_{i=1}^N b_i B_i(\vec r), \end{equation} where $b_i$ are the bending modes coefficients. The bending modes is an orthogonal set of vector, consequently the coefficients $b_i$ are given by \begin{equation} \label{eq:22} b_i = { \int \text{d}\vec r S(\vec r)B_i(\vec r) \over \int \text{d}\vec r \left| B_i(\vec r) \right|^2 }. \end{equation} \index{gmtMirrors!bending\_modes\_surface} <>= __device__ rtd bending_modes_surface(vector *v, const int N_mode, rtd *b, double *BM, const int BM_N_SAMPLE, const double BM_radius) { int ndx; rtd scale, S, s, t, fs, ft, onemt, onems; rtd *z; scale = 0.5/BM_radius; s = (scale*v->x + 0.5)*(BM_N_SAMPLE-1); t = (scale*v->y + 0.5)*(BM_N_SAMPLE-1); <> z = BM; S = ( ( z[ndx]*onemt + z[ndx+1]*t )*onems + ( z[ndx+BM_N_SAMPLE]*onemt + z[ndx+BM_N_SAMPLE+1]*t )*s ); return -S; } @ \index{gmtMirrors!partial\_x\_bending\_modes\_surface} <>= __device__ rtd partial_x_bending_modes_surface(vector *v, const int N_mode, rtd *b, double *BM, const int BM_N_SAMPLE, const double BM_radius) { int ndx, ndx_a, ndx_b; rtd scale, S, S_, s, t, fs, ft, onemt, onems, s_a, s_b, onems_a, onems_b, h; rtd *z; scale = 0.5/BM_radius; s = (scale*v->x + 0.5)*(BM_N_SAMPLE-1) + 1; t = (scale*v->y + 0.5)*(BM_N_SAMPLE-1); <> ndx_a = ndx; s_a = s; onems_a = onems; s = (scale*v->x + 0.5)*(BM_N_SAMPLE-1) - 1; t = (scale*v->y + 0.5)*(BM_N_SAMPLE-1); <> ndx_b = ndx; s_b = s; onems_b = onems; h = BM_radius*2/(BM_N_SAMPLE-1); z = BM; S_ = ( ( z[ndx_a]*onemt + z[ndx_a+1]*t )*onems_a + ( z[ndx_a+BM_N_SAMPLE]*onemt + z[ndx_a+BM_N_SAMPLE+1]*t )*s_a ); S_ -= ( ( z[ndx_b]*onemt + z[ndx_b+1]*t )*onems_b + ( z[ndx_b+BM_N_SAMPLE]*onemt + z[ndx_b+BM_N_SAMPLE+1]*t )*s_b ); S_ /= 2*h; S = S_; return -S; } @ \index{gmtMirrors!partial\_y\_bending\_modes\_surface} <>= __device__ rtd partial_y_bending_modes_surface(vector *v, const int N_mode, rtd *b, double *BM, const int BM_N_SAMPLE, const double BM_radius) { int ndx, ndx_a, ndx_b; rtd scale, S, S_, s, t, fs, ft, onemt, onems, t_a, t_b, onemt_a, onemt_b, h; rtd *z; scale = 0.5/BM_radius; s = (scale*v->x + 0.5)*(BM_N_SAMPLE-1); t = (scale*v->y + 0.5)*(BM_N_SAMPLE-1) + 1; <> ndx_a = ndx; t_a = t; onemt_a = onemt; s = (scale*v->x + 0.5)*(BM_N_SAMPLE-1); t = (scale*v->y + 0.5)*(BM_N_SAMPLE-1) - 1; <> ndx_b = ndx; t_b = t; onemt_b = onemt; h = BM_radius*2/(BM_N_SAMPLE-1); z = BM; S_ = ( ( z[ndx_a]*onemt_a + z[ndx_a+1]*t_a )*onems + ( z[ndx_a+BM_N_SAMPLE]*onemt_a + z[ndx_a+BM_N_SAMPLE+1]*t_a )*s ); S_ -= ( ( z[ndx_b]*onemt_b + z[ndx_b+1]*t_b )*onems + ( z[ndx_b+BM_N_SAMPLE]*onemt_b + z[ndx_b+BM_N_SAMPLE+1]*t_b )*s ); S_ /= 2*h; S = S_; return -S; } @ \index{gmtMirrors!partial\_y\_bending\_modes\_surface} <>= __device__ rtd partial_z_bending_modes_surface(void) { return 0.0; } @ \section{Modes} \label{sec:modes} The [[modes]] structure is a generalization of the the [[bending_modes]] structure; The [[bending_modes]] structure takes a single modal basis and applies to each segment. The [[modes]] structures has [[N_SET]] modal bases. The number of surface to be made of the bases set is [[N]]. If $[[N_SET]]=1$, then each surface uses exactly the same modal basis. But if $[[N_SET]]=[[N]]$, then each surface uses a different modal basis. \index{gmtMirrors!modes} <>= struct modes { <> <> void setup(char *filename, int N, int n_mode); void setup(int _BM_N_SAMPLE_, double _BM_radius_, int _N_SET_, int _N_MODE_, int *s2b, double *buffer, int _N_, int _n_mode_); void cleanup(void); void load(void); void load(int _BM_N_SAMPLE_, double _BM_radius_, int _N_SET_, int _N_MODE_, int *s2b, double *buffer); void reset_modes(double *buffer); void update(rtd *b); }; @ The parameters are \begin{itemize} \item the filename where the modes are saved: <>= char filename[256]; @ \item the number of modal basis set <>= int N_SET; @ \item the number of modes in the load modal basis, as opposed to [[n_mode]] the number of modes used to build a surface <>= int N_MODE; @ \item the index array of size [[N]] indicating for each surface which modal basis to use: <>= int *d__s2b; @ \end{itemize} @ \subsection{Setup \& Cleanup} \label{sec:setup--cleanup} The bending modes are loaded with: \index{gmtMirrors!modes!setup} <>= void modes::setup(char *_filename_, int _N_, int _n_mode_) { N = _N_; n_mode = _n_mode_; strcpy(filename,_filename_); load(); <> <> } void modes::setup(int _BM_N_SAMPLE_, double _BM_radius_, int _N_SET_, int _N_MODE_, int *s2b, double *buffer, int _N_, int _n_mode_) { N = _N_; n_mode = _n_mode_; load(_BM_N_SAMPLE_, _BM_radius_, _N_SET_, _N_MODE_, s2b, buffer); <> <> } @ Memory is freed with \index{gmtMirrors!modes!cleanup} <>= void modes::<> @ \subsection{Loading the modes} \label{sec:load-modes} The modes must be defined on a regular grid which size is restricted to a single segment. The linear grid sampling [[BM_N_SAMPLE]] and size [[BM_RADIUS]] is the same for all the segment. The number of modes per segment [[N_MODE]] is also the same. \subsubsection{From a binary file} \label{sec:bin-file} The binary file [[filename]] that contains the modes has the following format : [[N_SAMPLE]]:int$\times1$, [[LENGTH]]: double$\times1$, [[N_SET]]int$\times1$, [[N_MODE]]int$\times1$, [[d__s2b]]int$\times [[N]]$, [[d__BM]]: double$\times(7[[n_mode]][[BM_N_SAMPLE]]^2)$. The surface \#k with $k=[0,\dots,[[N]]-1]$ will use the basis in [[d__BM]] that starts at $[[d__s2b]][k]*[[N_MODE]]*[[N_SAMPLE]]^2$. \index{gmtMirrors!modes!load} <>= void modes::load(void) { FILE * pFile; long lSize; double * buffer; int *s2b; size_t result; //char path[256]; // sprintf(path, "%sgmtMirrors/%s.ceo",getenv("CEOPATH"),filename); INFO("@(CEO)>modes: loading modes from %s\n",filename); pFile = fopen ( filename, "rb" ); if (pFile==NULL) {fputs ("File error",stderr); exit (1);} // read modal basis info result = fread (&BM_N_SAMPLE,sizeof(int) ,1,pFile); result = fread (&BM_radius ,sizeof(double),1,pFile); result = fread (&N_SET,sizeof(int) ,1,pFile); result = fread (&N_MODE,sizeof(int) ,1,pFile); INFO("@(CEO)>modes: grid sampling and size: %d pixel, %f meter\n",BM_N_SAMPLE,BM_radius); INFO("@(CEO)>modes: number of set: %d and number of modes per set: %d\n",N_SET,N_MODE); // read surface 2 mode array s2b = (int*) malloc (N*sizeof(int)); result = fread (s2b ,sizeof(int),N,pFile); HANDLE_ERROR( cudaMalloc((void**)&d__s2b, N*sizeof(int)) ); HANDLE_ERROR( cudaMemcpy( d__s2b , s2b, N*sizeof(int), cudaMemcpyHostToDevice ) ); // allocate memory to contain the whole file: lSize = N_SET*N_MODE*BM_N_SAMPLE*BM_N_SAMPLE; buffer = (double*) malloc (lSize*sizeof(double)); if (buffer == NULL) {fputs ("Memory error",stderr); exit (2);} // copy the file into the buffer: result = fread (buffer ,sizeof(double),lSize,pFile); if (result != lSize) {fputs ("Reading error",stderr); exit (3);} /* the whole file is now loaded in the memory buffer. */ BM_radius *= 0.5; HANDLE_ERROR( cudaMalloc((void**)&d__BM_buffer, lSize*sizeof(double)) ); HANDLE_ERROR( cudaMemcpy( d__BM_buffer , buffer, lSize*sizeof(double), cudaMemcpyHostToDevice ) ); d__BM = d__BM_buffer; // terminate fclose (pFile); free(s2b); free (buffer); } @ \subsubsection{From a data set} \label{sec:data-set} The modes can also be loaded from input arrays: <>= void modes::load(int _BM_N_SAMPLE_, double _BM_radius_, int _N_SET_, int _N_MODE_, int *s2b, double *buffer) { long lSize; BM_N_SAMPLE = _BM_N_SAMPLE_; BM_radius = _BM_radius_; N_SET = _N_SET_; N_MODE = _N_MODE_; HANDLE_ERROR( cudaMalloc((void**)&d__s2b, N*sizeof(int)) ); HANDLE_ERROR( cudaMemcpy( d__s2b , s2b, N*sizeof(int), cudaMemcpyHostToDevice ) ); lSize = N_SET*N_MODE*BM_N_SAMPLE*BM_N_SAMPLE; BM_radius *= 0.5; HANDLE_ERROR( cudaMalloc((void**)&d__BM_buffer, lSize*sizeof(double)) ); HANDLE_ERROR( cudaMemcpy( d__BM_buffer , buffer, lSize*sizeof(double), cudaMemcpyHostToDevice ) ); d__BM = d__BM_buffer; } @ \subsection{Resetting the mode shapes} \label{sec:reset-shapes} \index{gmtMirrors!modes!reset\_modes} <>= void modes::reset_modes(double *buffer) { long lSize = N_SET*N_MODE*BM_N_SAMPLE*BM_N_SAMPLE; HANDLE_ERROR( cudaMemcpy( d__BM_buffer , buffer, lSize*sizeof(double), cudaMemcpyHostToDevice ) ); d__BM = d__BM_buffer; dim3 blockDim(N_THREAD,N_THREAD); dim3 gridDim(BM_N_SAMPLE/N_THREAD+1,BM_N_SAMPLE/N_THREAD+1,N); modes2surface LLL gridDim, blockDim RRR (d__BMS, n_mode, d__b, d__s2b, d__BM, N_MODE, BM_N_SAMPLE); } @ \subsection{Update} \label{sec:update} The bending mode coefficients [[d__b]] and the associated surface [[d__BMS]] are updated with: <>= void modes::update(rtd *b) { int n_byte = sizeof(rtd)*n_mode*N; HANDLE_ERROR( cudaMemcpy( d__b, b, n_byte, cudaMemcpyHostToDevice ) ); dim3 blockDim(N_THREAD,N_THREAD); dim3 gridDim(BM_N_SAMPLE/N_THREAD+1,BM_N_SAMPLE/N_THREAD+1,N); modes2surface LLL gridDim, blockDim RRR (d__BMS, n_mode, d__b, d__s2b, d__BM, N_MODE, BM_N_SAMPLE); } @ with the kernel <>= __global__ void modes2surface(double *BMS, const int n_mode, rtd *b, int *s2b, double *BM, const int N_MODE, const int N_SAMPLE) { int i, j, k, kz, l, k_mode, o; rtd *z; i = blockIdx.x * blockDim.x + threadIdx.x; j = blockIdx.y * blockDim.y + threadIdx.y; l = blockIdx.z; if ( (i>= struct stereoscopic_edge_sensors { <> void setup(gmt_m1 *_mirror_); void cleanup(void); void data(void); }; @ with \index{gmtMirrors!stereoscopic\_edge\_sensors!setup} \index{gmtMirrors!stereoscopic\_edge\_sensors!cleanup} \index{gmtMirrors!stereoscopic\_edge\_sensors!data} <>= void stereoscopic_edge_sensors::setup(gmt_m1 *_mirror_) { <> } void stereoscopic_edge_sensors::cleanup(void) { INFO("@(CEO)>stereoscopic_edge_sensors: freeing memory!\n"); <> } void stereoscopic_edge_sensors::data(void) { <> } @ The total number of unique edge sensors coordinates in the rigid body CS is 18 whose 12 are on the central segment and 6 on the peripheral segments leading to a total of 48 edge sensors coordinates in the GCS and 24 edge sensor mesurements. <>= int N, N_DATA; int n_byte, n_byte_DATA; <>= N = 48; N_DATA = 24; n_byte = sizeof(vector)*N; n_byte_DATA = sizeof(vector)*N_DATA; @ The $\alpha$ angle is set to 19.5degree <>= rtd alpha; <>= alpha = 19.5; @ The mirror the edge sensors are affected to is <>= gmt_m1 *mirror; <>= mirror = _mirror_; @ The coordinates of the edge sensors in the GCS are saved in the vectors [[d__v0]] for a perfectly aligned telescope and in the vectors [[d__v]] for a perturbed telescope. <>= vector *v0, *v, *d__v0, *d__v; <>= v0 = (vector *)malloc(n_byte); v = (vector *)malloc(n_byte); HANDLE_ERROR( cudaMalloc((void**)&d__v0, n_byte ) ); HANDLE_ERROR( cudaMalloc((void**)&d__v, n_byte ) ); <>= free( v0 ); free( v ); HANDLE_ERROR( cudaFree( d__v0 ) ); HANDLE_ERROR( cudaFree( d__v ) ); @ The edge sensor 24 measurements are saved in the vector [[d__dv0]] and [[d__dv]] for a perfectly aligned and perturbed telescope, respectively. <>= vector *dv0, *d__dv0; vector *dv, *d__dv; <>= dv0 = (vector *)malloc(n_byte_DATA); HANDLE_ERROR( cudaMalloc((void**)&d__dv0, n_byte_DATA ) ); dv = (vector *)malloc(n_byte_DATA); HANDLE_ERROR( cudaMalloc((void**)&d__dv, n_byte_DATA ) ); <>= free( dv0 ); HANDLE_ERROR( cudaFree( d__dv0 ) ); free( dv ); HANDLE_ERROR( cudaFree( d__dv ) ); @ The edge sensor coordinates in the GCS are computed with: <>= INFO("@(CEO)>gmt_m1: Initializing edge sensors coordinates!\n"); dim3 blockDim(6,1); dim3 gridDim(1,1); edge_sensors_read LLL gridDim, blockDim RRR (d__v0, mirror->N-1, mirror->D_full, alpha, mirror->rigid_body_CS.d__R, mirror->rigid_body_CS.d__origin, mirror->motion_CS.d__R, mirror->motion_CS.d__origin); HANDLE_ERROR( cudaMemcpy( v0, d__v0, n_byte, cudaMemcpyDeviceToHost ) ); HANDLE_ERROR( cudaMemcpy( d__v, d__v0, n_byte, cudaMemcpyDeviceToDevice ) ); HANDLE_ERROR( cudaMemcpy( v, d__v, n_byte, cudaMemcpyDeviceToHost ) ); @ and the reference measurements with: <>= HANDLE_ERROR( cudaMemset(d__dv0, 0, n_byte_DATA ) ); edge_sensors_data_kernel LLL gridDim, blockDim RRR (d__dv, d__v, d__dv0, mirror->N-1); HANDLE_ERROR( cudaMemcpy( d__dv0, d__dv, n_byte_DATA, cudaMemcpyDeviceToDevice ) ); HANDLE_ERROR( cudaMemcpy( dv0, d__dv0, n_byte_DATA, cudaMemcpyDeviceToHost ) ); HANDLE_ERROR( cudaMemset(d__dv, 0, n_byte_DATA ) ); HANDLE_ERROR( cudaMemcpy( dv, d__dv, n_byte_DATA, cudaMemcpyDeviceToHost ) ); @ with <>= __global__ void edge_sensors_read(vector *v, const int N_OUT_SEGMENT, const rtd D, const rtd alpha, rtd *d__rigid_body_R, vector *d__rigid_body_origin, rtd *d__motion_R, vector *d__motion_origin) { int i, iCoordSys, k; rtd R, a, offset_angle; i = threadIdx.x; if (i> __syncthreads(); <> } } @ First the coordinates $l_{7,k}$ are computed, <>= iCoordSys = 6; k = i; offset_angle = (3.0-2.0*i)/6.0 + a; <> @ followed by the coordinates $r_{7,k}$ <>= k += N_OUT_SEGMENT; offset_angle = (3.0-2.0*i)/6.0 - a; <> @ The coordinates of the edge sensors of segment \# $k=[[iCoordSys]]$ are computed next, starting with $l_{k,4}$ and $r_{k,4}$: <