8#ifndef __influenceFunctions_hpp__
9#define __influenceFunctions_hpp__
11#include "../improc/eigenCube.hpp"
12#include "../math/func/gaussian.hpp"
13#include "../ioutils/fits/fitsFile.hpp"
15#include <mx/math/eigenLapack.hpp>
17#include <mx/ioutils/pout.hpp>
18#include <mx/math/randomT.hpp>
30template <
typename realT>
31struct influenceFunctionGaussianSpec
41 std::vector<int> badActi;
42 std::vector<int> badActj;
43 std::vector<int> badActType;
45 influenceFunctionGaussianSpec()
56 realT pupilScale = pupilSz / dmSz;
58 realT act_space = ( dmSz - 1.0 ) / linNAct;
60 std::vector<realT> xcoords, ycoords;
62 realT xc, yc, xcen, ycen;
66 xcen = 0.5 * ( pupilSz - 1 );
67 ycen = 0.5 * ( pupilSz - 1 );
72 for(
int i = 0; i < linNAct; ++i )
74 xc = ( i + 0.5 ) * act_space - 0.5 * ( dmSz - 1 );
76 for(
int j = 0; j < linNAct; ++j )
80 if( badActi.size() > 0 && badActi.size() == badActj.size() )
82 for(
int k = 0;
k < badActi.size(); ++
k )
84 if( badActi[k] == i && badActj[k] == j )
94 yc = ( j + 0.5 ) * act_space - 0.5 * ( dmSz - 1 );
101 xoff = 0.5 * act_space;
104 _rad = sqrt( pow( xc + xoff, 2 ) + pow( yc, 2 ) );
106 if( _rad / act_space <= 0.5 * diameter + .01 || diameter == 0 )
108 xcoords.push_back( xcen + ( xc + xoff ) + 0.0 * ( dmSz - pupilSz ) );
109 ycoords.push_back( ycen + yc + 0.0 * ( dmSz - pupilSz ) );
115 std::cout <<
"Number of actuators: " << q <<
" " << xcoords.size() <<
"\n";
123 realT fwhm = act_space * 2. * sqrt( log2( 2 ) ) / sqrt( -log2( coupling ) );
127 for(
int i = 0; i < xcoords.size(); ++i )
133 act_inf.image( i ).data(), pupilSz, pupilSz, 0.0, 1.0, xc, yc,
fwhm2sigma( fwhm ) );
136 if( badActi.size() == badActj.size() && badActi.size() == badActType.size() )
138 Eigen::Array<realT, -1, -1> im;
140 std::cerr <<
"ok trying bad acts\n";
141 std::cerr << badActi.size() <<
"\n";
142 for(
int k = 0;
k < badActi.size(); ++
k )
144 if( badActType[k] == 0 )
147 im.resize( pupilSz, pupilSz );
149 xc = ( badActi[
k] + 0.5 ) * act_space - 0.5 * ( dmSz - 1 );
150 yc = ( badActj[
k] + 0.5 ) * act_space - 0.5 * ( dmSz - 1 );
157 xoff = 0.5 * act_space;
160 _rad = sqrt( pow( xc + xoff, 2 ) + pow( yc, 2 ) );
162 std::cerr << xc <<
" " << yc <<
"\n";
164 if( _rad / act_space <= 0.5 * diameter + .01 || diameter == 0 )
166 xc = ( xcen + ( xc + xoff ) + 0.0 * ( dmSz - pupilSz ) );
167 yc = ( ycen + yc + 0.0 * ( dmSz - pupilSz ) );
169 std::cerr << xc <<
" " << yc <<
"\n";
173 for(
int i = 0; i < xcoords.size(); ++i )
175 act_inf.image( i ) += im;
182 std::cerr <<
"Didn't do bad atuators!\n";
188 head.
append(
"LINNACT", linNAct,
"Linear number of actuators" );
189 head.
append(
"DIAMACT", diameter,
"Diameter of DM in actuators" );
190 head.
append(
"DIAMPIX", diameter * act_space,
"Diameter of DM in pixels" );
191 head.
append(
"ACTSPACE", act_space,
"Actuator spacing, in pixels" );
192 head.
append(
"IFTYPE",
"GAUSSIAN",
"Type of influence function" );
193 head.
append(
"COUPLING", coupling,
"Coupling at 1 actuator spacing" );
194 head.
append(
"FWHM", fwhm / act_space,
"Resultant FWHM in actuators" );
195 head.
append(
"OFFSTODD", offsetOdd,
"Whether odd rows are offset by 1/2" );
198 ff.
write( fName, act_inf, head );
203 Eigen::Array<realT, -1, -1> coords( 2, xcoords.size() );
204 for(
int i = 0; i < xcoords.size(); ++i )
206 coords( 0, i ) = xcoords[i];
207 coords( 1, i ) = ycoords[i];
211 ff.
write( fName, coords );
229template <
typename realT>
231 const std::string &dmName,
238 bool offsetOdd =
false )
243 realT pupilScale = pupilSz / dmSz;
245 realT act_space = ( dmSz - 1.0 ) / linNAct;
247 std::vector<realT> xcoords, ycoords;
249 realT xc, yc, xcen, ycen;
253 xcen = 0.5 * ( pupilSz - 1 );
254 ycen = 0.5 * ( pupilSz - 1 );
260 if( couplingRange > 0 )
265 for(
int i = 0; i < linNAct; ++i )
267 xc = ( i + 0.5 ) * act_space - 0.5 * ( dmSz - 1 );
269 for(
int j = 0; j < linNAct; ++j )
271 yc = ( j + 0.5 ) * act_space - 0.5 * ( dmSz - 1 );
278 xoff = 0.5 * act_space;
281 _rad = sqrt( pow( xc + xoff, 2 ) + pow( yc, 2 ) );
283 if( _rad / act_space <= 0.5 * diameter + .01 || diameter == 0 )
285 xcoords.push_back( xcen + ( xc + xoff ) + 0.0 * ( dmSz - pupilSz ) );
286 ycoords.push_back( ycen + yc + 0.0 * ( dmSz - pupilSz ) );
291 std::cout <<
"Number of actuators: " << xcoords.size() <<
"\n";
300 realT fwhm = act_space * 2. * sqrt( log2( 2 ) ) / sqrt( -log2( coupling ) );
306 for(
int i = 0; i < xcoords.size(); ++i )
311 if( couplingRange > 0 )
313 fw = fwhm + fwhm * ( 1.0 - 2.0 * uniVar ) * couplingRange;
320 mx::math::func::gaussian2D<realT>(
327 head.
append(
"LINNACT", linNAct,
"Linear number of actuators" );
328 head.
append(
"DIAMACT", diameter,
"Diameter of DM in actuators" );
329 head.
append(
"DIAMPIX", diameter * act_space,
"Diameter of DM in pixels" );
330 head.
append(
"ACTSPACE", act_space,
"Actuator spacing, in pixels" );
331 head.
append(
"IFTYPE",
"GAUSSIAN",
"Type of influence function" );
332 head.
append(
"COUPLING", coupling,
"Coupling at 1 actuator spacing" );
333 head.
append(
"FWHM", fwhm / act_space,
"Resultant FWHM in actuators" );
334 head.
append(
"OFFSTODD", offsetOdd,
"Whether odd rows are offset by 1/2" );
336 std::string fName = mx::AO::path::dm::influenceFunctions( dmName,
true );
337 ff.
write( fName, act_inf, head );
342 Eigen::Array<realT, -1, -1> coords( 2, xcoords.size() );
343 for(
int i = 0; i < xcoords.size(); ++i )
345 coords( 0, i ) = xcoords[i];
346 coords( 1, i ) = ycoords[i];
349 fName = mx::AO::path::dm::actuatorPositions( dmName,
true );
350 ff.
write( fName, coords );
365template <
typename realT>
366void ifPInv(
const std::string &dmName, realT maxCondition = -1 )
368 std::string pinvName = mx::AO::path::dm::pseudoInverse( dmName );
369 std::string mmodesName = mx::AO::path::dm::mirrorModes( dmName,
true );
370 std::string svsName = mx::AO::path::dm::singularValues( dmName );
371 std::string ifName = mx::AO::path::dm::influenceFunctions( dmName );
373 std::cout <<
"mmodesName: " << mmodesName <<
"\n";
378 ff.
read( ifName, actInf );
380 int nrows = actInf.rows();
381 int ncols = actInf.cols();
382 int nplanes = actInf.planes();
384 Eigen::Array<realT, -1, -1> rowInf = actInf.
asVectors();
385 actInf.resize( 0, 0, 0 );
387 Eigen::Array<realT, -1, -1> PInv;
392 Eigen::Array<realT, -1, -1> U, S, VT;
394 int interact = MX_PINV_PLOT | MX_PINV_ASK;
395 if( maxCondition >= 0 )
396 interact = MX_PINV_NO_INTERACT;
401 head.
append(
"MAXCONDN", maxCondition,
"Max. condition no. in pseudo inverse" );
402 head.
append(
"CONDN", condition,
"Actual condition number of pseudo inverse" );
403 head.
append(
"NREJECT", nRejected,
"Number of s.v.s rejected" );
407 ff.
write( pinvName, PInv, head );
412 ff.
write( mmodesName, mmodes, head );
415 fout.open( svsName );
417 for(
int i = 0; i < S.rows(); ++i )
419 fout << S( i ) <<
"\n";
366void ifPInv(
const std::string &dmName, realT maxCondition = -1 ) {
…}
452template <
typename realT>
455 M2c = Ap.matrix().transpose() * M.
asVectors().matrix();
469template <
typename realT>
470void m2cMatrix(
const std::string &dmName,
const std::string &basisName )
472 std::string M2cFName;
475 Eigen::Array<realT, -1, -1> Ap;
476 std::string pinvFName = mx::AO::path::dm::pseudoInverse( dmName );
477 ff.
read( pinvFName, Ap );
480 std::string basisFName;
482 basisFName = mx::AO::path::basis::modes( basisName );
483 ff.
read( basisFName, M );
485 Eigen::Array<realT, -1, -1> M2c;
488 M2cFName = mx::AO::path::dm::M2c( dmName, basisName,
true );
490 ff.
write( M2cFName, M2c );
470void m2cMatrix(
const std::string &dmName,
const std::string &basisName ) {
…}
501template <
typename realT>
508 std::string basisFName;
510 basisFName = mx::AO::path::basis::modes( basisName );
512 ff.
read( basisFName, M );
514 Eigen::Array<realT, -1, -1> M2c;
516 M2c.resize( M.planes(), M.planes() );
519 for(
int i = 0; i < M2c.rows(); ++i )
524 std::string M2cFName;
526 M2cFName = mx::AO::path::dm::M2c(
"modalDM", basisName,
true );
528 ff.
write( M2cFName, M2c );
539template <
typename realT>
540void modesOnDM(
const std::string &dmName,
const std::string &basisName )
542 Eigen::Array<realT, -1, -1> m2c;
546 m2cName = mx::AO::path::dm::M2c( dmName, basisName );
549 ff.
read( m2cName, m2c );
551 std::cout << m2c.rows() <<
"\n";
552 std::cout << m2c.cols() <<
"\n";
554 int nmodes = m2c.cols();
555 int nact = m2c.rows();
561 std::string fName = mx::AO::path::dm::influenceFunctions( dmName,
true );
562 ff.
read( fName, act_inf );
564 projModes.resize( act_inf.rows(), act_inf.cols(), nmodes );
567#pragma omp parallel for
568 for(
int i = 0; i < nmodes; ++i )
570 Eigen::Array<realT, -1, -1> m, c;
572 std::cout << i + 1 <<
"/" << nmodes <<
"\n";
573 m.resize( nmodes, 1 );
577 c = m2c.matrix() * m.matrix();
585 for(
int j = 0; j < nact; ++j )
587 projModes.
image( i ) += c( j, 0 ) * act_inf.
image( j );
596 oName = mx::AO::path::dm::projectedModes( dmName, basisName,
true );
598 ff.
write( oName, projModes );
540void modesOnDM(
const std::string &dmName,
const std::string &basisName ) {
…}
Standardized paths for the mx::AO system.
std::string actuatorPositions(const std::string &dmName, bool create=false)
The path for the deformable mirror (DM) actuator positions.
std::string influenceFunctions(const std::string &dmName, bool create=false)
The path for the deformable mirror (DM) influence functions.
Class to manage interactions with a FITS file.
int write(const dataT *im, int d1, int d2, int d3, fitsHeader *head)
Write the contents of a raw array to the FITS file.
int read(dataT *data)
Read the contents of the FITS file into an array.
An image cube with an Eigen-like API.
Eigen::Map< Eigen::Array< dataT, Eigen::Dynamic, Eigen::Dynamic > > asVectors()
Return an Eigen::Eigen::Map of the cube where each image is a vector.
Eigen::Map< Eigen::Array< dataT, Eigen::Dynamic, Eigen::Dynamic > > image(Index n)
Returns a 2D Eigen::Eigen::Map pointed at the specified image.
A random number type, which functions like any other arithmetic type.
void seed(typename ranengT::result_type seedval)
Set the seed of the random engine.
constexpr units::realT k()
Boltzmann Constant.
int eigenPseudoInverse(Eigen::Array< dataT, -1, -1 > &PInv, dataT &condition, int &nRejected, Eigen::Array< dataT, -1, -1 > &U, Eigen::Array< dataT, -1, -1 > &S, Eigen::Array< dataT, -1, -1 > &VT, Eigen::Array< dataT, -1, -1 > &A, dataT &maxCondition, dataT alpha=0, int interact=MX_PINV_NO_INTERACT)
Calculate the pseudo-inverse of a patrix using the SVD.
realT gaussian2D(const realT x, const realT y, const realT G0, const realT G, const realT x0, const realT y0, const realT sigma)
Find value at position (x,y) of the 2D arbitrarily-centered symmetric Gaussian.
floatT fwhm2sigma(floatT fw)
Convert from FWHM to the Gaussian width parameter.
void modalDMM2cMatrix(const std::string &basisName)
Generate the modes-to-commands matrix for a set of modes on the modal DM.
void modesOnDM(const std::string &dmName, const std::string &basisName)
Calculate the basis set as it will be reproduced by the DM.
void m2cMatrix(Eigen::Array< realT, -1, -1 > &M2c, Eigen::Array< realT, -1, -1 > &Ap, mx::improc::eigenCube< realT > &M)
Calculate the modes-to-commands matrix for a set of modes.
void ifPInv(const std::string &dmName, realT maxCondition=-1)
Calculate the pseudo-inverse and mirror modes for a set of influence functions.
void influenceFunctionsGaussian(const std::string &dmName, realT dmSz, int linNAct, realT diameter, realT coupling, realT couplingRange, realT pupilSz=0, bool offsetOdd=false)
Create a set of Gaussian influence functions (IFs) on a square grid.