fourierCovariance.hpp

Go to the documentation of this file.

/** \file fourierCovariance.hpp
  * \author Jared R. Males (jaredmales@gmail.com)
  * \brief Calculation of the modal covariance in the Fourier basis.
  * \ingroup mxAO_files
  *
  */
 
#ifndef fourierCovariance_hpp
#define fourierCovariance_hpp
 
 
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
 
#include "../../math/constants.hpp"
#include "../../math/func/jinc.hpp"
#include "../../sigproc/fourierModes.hpp"
#include "../../ioutils/fits/fitsFile.hpp"
#include "../../improc/eigenImage.hpp"
#include "../../improc/eigenCube.hpp"
#include "../../mxlib_uncomp_version.h"
#include "../../ipc/ompLoopWatcher.hpp"
#include "../../sys/timeUtils.hpp"
#include "../../math/eigenLapack.hpp"
 
#include "../../math/func/airyPattern.hpp"
 
 
 
#include "aoAtmosphere.hpp"
#include "aoPSDs.hpp"
#include "aoSystem.hpp"
#include "varmapToImage.hpp"
 
 
namespace mx
{
namespace AO
{
namespace analysis
{
 
#ifndef WSZ
 
/** \def WSZ
  * \brief The size of the gsl_integration workspaces.
  */
#define WSZ 100000
 
#endif
 
 
//Forward Declaration
template<typename realT, typename aosysT>
struct fourierCovariance;
 
///Worker for the azimuthal integral (in phi) for the basic Fourier mode covariance.
/**
  * \param phi the angle at which to evaluate the integrand
  * \param params a pointer to a object of type fourierCovariance<realT, aosyT>
  *
  * \tparam realT a floating point type used for all calculations
  * \tparam aosysT the type of the AO system structure
  *
  * \ingroup mxAOAnalytic
  */
template<typename realT, typename aosysT>
realT phiInt_basic (realT phi, void * params)
{
 
   fourierCovariance<realT, aosysT> * Pp = (fourierCovariance<realT, aosysT> *) params;
 
   realT D = Pp->aosys->D();
 
   int p = Pp->p;
   realT m = Pp->m;
   realT n = Pp->n;
 
   int pp = Pp->pp;
   realT mp = Pp->mp;
   realT np = Pp->np;
 
   realT k = Pp->k;
 
   realT cosp = cos(phi);
   realT sinp = sin(phi);
 
   /*** no prime ***/
   realT kmn_p, kmn_m;
   realT Ji_mn_p, Ji_mn_m;
   realT Q_mn;//, Qs_mn;
 
   kmn_p = sqrt( pow(k*cosp + m/D, 2) + pow(k*sinp + n/D, 2));
   kmn_m = sqrt( pow(k*cosp - m/D, 2) + pow(k*sinp - n/D, 2));
 
   Ji_mn_p = math::func::jinc(math::pi<realT>()*D*kmn_p);
   Ji_mn_m = math::func::jinc(math::pi<realT>()*D*kmn_m);
 
 
   realT N = 1./sqrt(0.5 + p*math::func::jinc(math::two_pi<realT>()*sqrt(m*m+n*n)));
 
   Q_mn = N*(Ji_mn_p + p*Ji_mn_m);
 
   /*** primed ***/
   realT kmpnp_p, kmpnp_m;
   realT Ji_mpnp_p, Ji_mpnp_m;
   realT Q_mpnp;//, Qs_mpnp;
 
   kmpnp_p = sqrt( pow(k*cosp + mp/D, 2) + pow(k*sinp + np/D, 2));
   kmpnp_m = sqrt( pow(k*cosp - mp/D, 2) + pow(k*sinp - np/D, 2));
 
 
   Ji_mpnp_p = math::func::jinc(math::pi<realT>()*D*kmpnp_p);
   Ji_mpnp_m = math::func::jinc(math::pi<realT>()*D*kmpnp_m);
 
   realT Np = 1./sqrt(0.5 + pp*math::func::jinc(math::two_pi<realT>()*sqrt(mp*mp+np*np)));
 
   Q_mpnp = Np*(Ji_mpnp_p + pp*Ji_mpnp_m);
 
   realT P = Pp->aosys->psd(Pp->aosys->atm, 0, k, 1.0); //vonKarmanPSD(k, D, L0, Pp->subPiston, Pp->subTipTilt);
 
   return P*k*(Q_mn*Q_mpnp);
}
 
///Worker for the azimuthal integral (in phi) for the modified Fourier mode covariance.
/**
  * \param phi the angle at which to evaluate the integrand
  * \param params a pointer to a object of type fourierCovariance<realT, aosyT>
  *
  * \tparam realT a floating point type used for all calculations
  * \tparam aosysT the type of the AO system structure
  */
template<typename realT, typename aosysT>
realT phiInt_mod (realT phi, void * params)
{
 
   fourierCovariance<realT, aosysT> * Pp = (fourierCovariance<realT, aosysT> *) params;
 
   //realT L0 = Pp->aosys->atm.L_0();
   realT D = Pp->aosys->D();
 
   int p = Pp->p;
   realT m = Pp->m;
   realT n = Pp->n;
 
   int pp = Pp->pp;
   realT mp = Pp->mp;
   realT np = Pp->np;
 
   realT k = Pp->k;
 
   realT mnCon = Pp->mnCon;
 
   realT cosp = cos(phi);
   realT sinp = sin(phi);
 
   /*** no prime ***/
   realT kmn_p, kmn_m;
   realT Ji_mn_p, Ji_mn_m;
   
   kmn_p = sqrt( pow(k*cosp + m/D, 2) + pow(k*sinp + n/D, 2));
   kmn_m = sqrt( pow(k*cosp - m/D, 2) + pow(k*sinp - n/D, 2));
 
   Ji_mn_p = math::func::jinc(math::pi<realT>()*D*kmn_p);
   Ji_mn_m = math::func::jinc(math::pi<realT>()*D*kmn_m);
 
 
   /*** primed ***/
   realT kmpnp_p, kmpnp_m;
   realT Ji_mpnp_p, Ji_mpnp_m;
   
   kmpnp_p = sqrt( pow(k*cosp + mp/D, 2) + pow(k*sinp + np/D, 2));
   kmpnp_m = sqrt( pow(k*cosp - mp/D, 2) + pow(k*sinp - np/D, 2));
 
 
   Ji_mpnp_p = math::func::jinc(math::pi<realT>()*D*kmpnp_p);
   Ji_mpnp_m = math::func::jinc(math::pi<realT>()*D*kmpnp_m);
 
   realT QQ;
 
   if(p == pp)
   {
      QQ = 2.0*( Ji_mn_p*Ji_mpnp_p + Ji_mn_m*Ji_mpnp_m);
   }
   else
   {
      QQ = 2.0*( Ji_mn_p*Ji_mpnp_m + Ji_mn_m*Ji_mpnp_p);
   }
 
   realT P = Pp->aosys->psd(Pp->aosys->atm, 0, k, 1.0); //psd(Pp->aosys->atm, k, Pp->aosys->lam_sci(), Pp->aosys->lam_wfs(), 1.0);
 
   /*if(mnCon > 0 )
   {
      if( k*D <= mnCon )
      {
         P *= pow(math::two_pi<realT>()*Pp->aosys->atm.v_wind()* k * (Pp->aosys->minTauWFS(0)+Pp->aosys->deltaTau()),2);
      }
   }*/
 
 
   return P*k*QQ;
}
 
 
 
///Worker function for the radial integral in the covariance calculation
/**
  * \param k the spatial frequency at which to evaluate the integrand
  * \param params a pointer to a object of type fourierCovariance<realT, aosyT>
  *
  * \tparam realT a floating point type used for all calculations.  As of Nov 2016 must be double due to gsl_integration.
  * \tparam aosysT the type of the AO system structure
  */
template<typename realT, typename aosysT>
realT kInt (realT k, void * params)
{
   fourierCovariance<realT, aosysT> * Pp = (fourierCovariance<realT, aosysT> *) params;
 
   realT result, error;
 
   gsl_function func;
 
   //Here choose between basic and modified Fourier modes.
   if(Pp->useBasic)
   {
      func.function = &phiInt_basic<realT, aosysT>;
   }
   else
   {
      func.function = &phiInt_mod<realT, aosysT>;
   }
   func.params = Pp;
 
   Pp->k = k;
 
   //Tolerances of 1e-4 seem to be all we can ask for
   gsl_integration_qag(&func, 0., math::two_pi<double>(), 1e-10, 1e-4, WSZ, GSL_INTEG_GAUSS21, Pp->phi_w, &result, &error);
 
 
   return result;
}
 
/// Structure to manage the Fourier mode covariance calculation, passed to integration functions
/**
  * \tparam realT a floating point type used for all calculations.  As of Nov 2016 must be double due to gsl_integration.
  * \tparam aosysT the type of the AO system structure
  */
template<typename realT, typename aosysT>
struct fourierCovariance
{
   ///Pointer to an AO system, which contains the relevant spatial PSD of turbulence.
   aosysT * aosys {nullptr};
 
   ///Flag controlling use of basic or modified Fourier modes.  If true, the basic sin/cos modes are used.  If false (default), the modified modes are u sed.
   bool useBasic {false};
 
   ///p-index of the unprimed mode.  +/-1 for modified modes.  If basic, then +1==>cosine, -1==>sine.
   int p;
 
   ///The m-index of the unprimed mode, corresponding to the \f$ k_u = m/D \f$ component of spatial frequency.
   realT m;
 
   ///The n-indexof the unprimed mode, corresponding to the \f$ k_v = n/D \f$ component of spatial frequency.
   realT n;
 
   ///p-index of the primed mode.  +/-1 for modified modes.  If basic, then +1==>cosine, -1==>sine.
   int pp;
 
   ///The m-index of the primed mode, corresponding to the \f$ k_u = m/D \f$ component of spatial frequency.
   realT mp;
 
   ///The n-indexof the primed mode, corresponding to the \f$ k_v = n/D \f$ component of spatial frequency.
   realT np;
 
   ///Spatial frequency being calculated, passed for use in the integrand worker functions.
   realT k;
 
   ///The maximum controlled value of spatial frequency (k*D).  If < 0 then not controlled.
   realT mnCon {0};
 
   ///Absolute tolerance for the radial integral.  Default is 1e-7.
   realT absTol {1e-7};
 
   ///Relative tolerance for the radial integral. Default is 1e-7.
   realT relTol {1e-7};
 
   ///Working memory for the azimuthal integral.
   gsl_integration_workspace * phi_w;
 
   ///Working memory for the radial integral.
   gsl_integration_workspace * k_w;
 
   ///Constructor
   fourierCovariance()
   {
      phi_w = gsl_integration_workspace_alloc (WSZ);
      k_w = gsl_integration_workspace_alloc (WSZ);
   }
 
   ///Destructor
   ~fourierCovariance()
   {
      gsl_integration_workspace_free(phi_w);
      gsl_integration_workspace_free(k_w);
   }
 
   ///Calculate the covariance between the two modes.
   /** \todo document me
     * \todo handle gsl errors
     */
   realT getVariance(realT & error)
   {
      realT result;
 
      gsl_function func;
      func.function = &kInt<realT, aosysT>;
      func.params = this;
 
      gsl_set_error_handler_off();
 
      int ec = gsl_integration_qagiu (&func, 0, absTol, relTol, WSZ, k_w, &result, &error);
 
      //Finally, convert to square sampling.
      return result;
   }
 
};
 
///Calculate a vector of Fourier mode variances.
template<typename realT, typename aosysT>
int fourierVarVec( const std::string & fname,
                   int N,
                   aosysT & aosys,
                   realT absTol,
                   realT relTol,
                   bool modifed=true
                 )
{
 
   std::vector<realT> var(N,0);
 
   ipc::ompLoopWatcher<> watcher(N, std::cerr);
 
   realT mnCon = 0;
   if( aosys.d_min(0) > 0)
   {
      mnCon = floor( aosys.D()/aosys.d_min(0)/2.0);
   }
 
   #pragma omp parallel
   {
      fourierCovariance<realT, aosysT> Pp;
      Pp.absTol = absTol;
      Pp.relTol = relTol;
      Pp.aosys = &aosys;
 
      Pp.mnCon = mnCon;
 
      realT result, error;
 
      #pragma omp for schedule(dynamic,5)
      for(int i=0; i< N; ++i)
      {
         Pp.p = +1;
         Pp.m = i+1;
         Pp.n = 0;
 
         Pp.pp = +1;
         Pp.mp = i+1;
         Pp.np = 0;
 
         result = Pp.getVariance(error);
 
         var[i] = result;
 
         watcher.incrementAndOutputStatus();
      }
   }
 
   std::ofstream fout;
   fout.open(fname);
 
   realT D = aosys.D();
   realT r_0 = aosys.atm.r_0();
   realT tot = 1.0299*pow(D/r_0, 5./3.);
 
   realT sum = 0;
 
   std::cout << 0 << " " << 0 << " " << 0 << " " <<  tot << "\n";
   for(int i=0; i<N; ++i)
   {
      realT k = (i+1)/D;
 
      realT P = aosys.psd(aosys.atm, 0, k, 1.0);// aosys.lam_sci(), 0, aosys.lam_wfs(), 1.0);
 
      if(mnCon > 0 )
      {
         if( k*D < mnCon )
         {
            P *= pow(math::two_pi<realT>()*aosys.atm.v_wind()* k * (aosys.minTauWFS(0)+aosys.deltaTau()),2);
         }
      }
 
      sum += var[i];
      fout << i+1 << " " << var[i] << " " << P/pow(D,2) << " " << tot-sum << "\n";
   }
 
   fout.close();
 
   return 0;
}
 
///Calculate a map of Fourier variances by convolution with the PSD
/** Uses the Airy pattern for the circularly unobstructed aperture.
  *
  * \returns 0 on success
  * \returns -1 on error
  */
template<typename realT, typename aosysT>
int fourierPSDMap( improc::eigenImage<realT> & var, ///< [out] The variance estimated by convolution with the PSD
                   improc::eigenImage<realT> & psd, ///< [out] the PSD map
                   int N, ///< [in] the number of components to analyze
                   int overSample,
                   aosysT & aosys ///< [in[ the AO system defining the PSD characteristics.
                 )
{
   N *= overSample;
   psd.resize(2*N + 1, 2*N+1);
 
   realT mnCon = 0;
   if( aosys.d_min() > 0)
   {
      mnCon = floor( aosys.D()/aosys.d_min()/2.0);
   }
 
   for(int i=0; i<=N; ++i)
   {
      for(int j=-N; j<=N; ++j)
      {
 
         realT D = aosys.D();
         realT k = sqrt( pow(i,2) + pow(j,2))/D/overSample;
 
         realT P = aosys.psd(aosys.atm, k, aosys.lam_sci(), 0, aosys.lam_wfs(), 1.0);
 
         if(mnCon > 0 )
         {
            if( k*D < mnCon )
            {
               P *= pow(math::two_pi<realT>()*aosys.atm.v_wind()* k * (aosys.minTauWFS()+aosys.deltaTau()),2);
            }
         }
 
         psd(N+i, N + j) = P/pow(D*overSample,2);
         psd(N-i, N - j) = P/pow(D*overSample,2);
      }
   }
 
   //Create Airy PSF for convolution with PSD psd.
   Eigen::Array<realT, -1,-1> psf;
   psf.resize(2*N+3,2*N+3);
   for(int i=0;i<psf.rows();++i)
   {
      for(int j=0;j<psf.cols();++j)
      {
         psf(i,j) = mx::math::func::airyPattern(sqrt( pow( i-floor(.5*psf.rows()),2) + pow(j-floor(.5*psf.cols()),2))/overSample);
      }
   }
 
   mx::AO::analysis::varmapToImage(var, psd, psf);
 
   return 0;
}
 
template<typename realT>
int fourierCovarMap( const std::string & fname, ///< [out] the path where the output FITS file will be written
                     int N, ///< [in] the linear number of Fourier modes across the aperture.  The Nyquist frequency is set by N/2.
                     realT D,
                     realT L_0,
                     bool subPist,
                     bool subTilt,
                     realT absTol,
                     realT relTol,
                     bool modified=true
                   )
{
   std::vector<mx::sigproc::fourierModeDef> ml;
   mx::sigproc::makeFourierModeFreqs_Rect(ml, N);
 
 
   aoSystem<realT, vonKarmanSpectrum<realT>> aosys;
 
   aosys.loadMagAOX();
 
   //This is just a normalization parameter in this context.
   aosys.atm.r_0(1.0, 0.5e-6);
 
   aosys.D( D );
   aosys.atm.L_0( L_0 );
   aosys.psd.subPiston( subPist );
   aosys.psd.subTipTilt( subTilt );
 
   int psz = ml.size();
 
   Eigen::Array<realT,-1,-1> covar( psz, psz);
   covar.setZero();
 
   //int ncalc = 0.5*( psz*psz - psz);
 
   ipc::ompLoopWatcher<> watcher(psz, std::cout);
 
   std::cerr << "Starting . . .\n";
   #pragma omp parallel
   {
      fourierCovariance<realT, aoSystem<realT, vonKarmanSpectrum<realT> > > Pp;
      Pp.absTol = absTol;
      Pp.relTol = relTol;
      Pp.aosys = &aosys;
 
      if(!modified) Pp.useBasic = true;
 
      realT result, error;
 
      #pragma omp for schedule(static,5)
      for(int i=0; i< psz; ++i)
      {
         for(int j=i; j< psz; ++j)
         {
            Pp.p = ml[i].p;
            Pp.m = ml[i].m;
            Pp.n = ml[i].n;
 
            Pp.pp = ml[j].p;
            Pp.mp = ml[j].m;
            Pp.np = ml[j].n;
            result = Pp.getVariance(error);
 
            covar(i,j) = result;
 
         }
         watcher.incrementAndOutputStatus();
      }
   }
 
   fits::fitsHeader head;
   head.append("DIAMETER", aosys.D(), "Diameter in meters");
   head.append("NSUBAP", N, "Linear number of s.f. sampled");
   head.append("L0", aosys.atm.L_0(0), "Outer scale (L_0) in meters");
   head.append("SUBPIST", aosys.psd.subPiston(), "Piston subtractioon true/false flag");
   head.append("SUBTILT", aosys.psd.subTipTilt(), "Tip/Tilt subtractioon true/false flag");
   head.append("ABSTOL", absTol, "Absolute tolerance in qagiu");
   head.append("RELTOL", relTol, "Relative tolerance in qagiu");
 
   fitsHeaderGitStatus(head, "mxlib_uncomp",  MXLIB_UNCOMP_CURRENT_SHA1, MXLIB_UNCOMP_REPO_MODIFIED);
   
 
   fits::fitsFile<realT> ff;
   ff.write(fname + ".fits", covar, head);
 
   return 0;
}
 
 
 
 
template<typename realT, typename aosysT>
int fourierCovarMapSeparated( const std::string & fname,
                              int N,
                              aosysT & aosys,
                              realT absTol,
                              realT relTol,
                              bool modified=true)
{
   std::vector<mx::sigproc::fourierModeDef> ml;
   mx::sigproc::makeFourierModeFreqs_Rect(ml, N);
 
   int psz = 0.5*ml.size();
 
   Eigen::Array<realT,-1,-1> covar_pp( (int) (0.5*psz), (int)(.5*psz)), covar_ppp( (int) (0.5*psz), (int)(0.5*psz));
   covar_pp.setZero();
   covar_ppp.setZero();
 
   realT mnCon = 0;
   if( aosys.d_min() > 0)
   {
      mnCon = floor( aosys.D()/aosys.d_min()/2.0);
   }
 
   ipc::ompLoopWatcher<> watcher((psz+1)*0.125*(psz+1)*2, std::cout);
 
   #pragma omp parallel
   {
 
      fourierCovariance<realT, aosysT > Pp;
      Pp.absTol = absTol;
      Pp.relTol = relTol;
      Pp.aosys = &aosys;
 
      if(!modified) Pp.useBasic = true;
 
      Pp.mnCon = mnCon;
 
      realT result, error;
 
      #pragma omp for schedule(dynamic,5)
      for(int i=0; i< psz; i+=2)
      {
         for(int j=0; j<= 0.5*i; ++j)
         {
            for(int k=0; k< 2; ++k)
            {
               Pp.p = ml[i].p;
               Pp.m = ml[i].m;
               Pp.n = ml[i].n;
 
               Pp.pp = ml[2*j + k].p;
               Pp.mp = ml[2*j + k].m;
               Pp.np = ml[2*j + k].n;
               result = Pp.getVariance(error);
 
               if( Pp.p == Pp.pp)
               {
                  covar_pp(i/2, j) = result;
               }
               else
               {
                  covar_ppp(i/2, j) = result;
               }
               watcher.incrementAndOutputStatus();
            }
         }
      }
   }
 
   fits::fitsHeader head;
   head.append("DIAMETER", aosys.D(), "Diameter in meters");
   head.append("L0", aosys.atm.L_0(), "Outer scale (L_0) in meters");
   head.append("SUBPIST", aosys.psd.subPiston(), "Piston subtractioon true/false flag");
   head.append("SUBTILT", aosys.psd.subTipTilt(), "Tip/Tilt subtractioon true/false flag");
   head.append("ABSTOL", absTol, "Absolute tolerance in qagiu");
   head.append("RELTOL", relTol, "Relative tolerance in qagiu");
 
   
   fitsHeaderGitStatus(head, "mxlib_uncomp",  MXLIB_UNCOMP_CURRENT_SHA1, MXLIB_UNCOMP_REPO_MODIFIED);
   
   fits::fitsFile<realT> ff;
   ff.write(fname + "_pp.fits", covar_pp, head);
   ff.write(fname + "_ppp.fits", covar_ppp, head);
 
   return 0;
}
 
template<typename realT>
void calcKLCoeffs( const std::string & outFile,
                   const std::string & cvFile )
{
   fits::fitsFile<realT> ff;
 
   Eigen::Array<realT,-1,-1> cvT, cv, evecs, evals;
 
   ff.read(cv, cvFile);
 
   //cvT = cv.block(0,0, 1000,1000);//.transpose();
 
   std::cerr << cvT.rows() << " " << cvT.cols() << "\n";
   std::cerr << "1\n";
   math::syevrMem<double> mem;
 
   double t0 = sys::get_curr_time();
 
   int info = math::eigenSYEVR<double,double>(evecs, evals, cv, 0, -1, 'U', &mem);
 
   double t1 = sys::get_curr_time();
 
   std::cerr << "2\n";
 
   if(info !=0 )
   {
      std::cerr << "info =" << info << "\n";
      exit(0);
   }
 
   std::cerr << "Time = " << t1-t0 << " secs\n";
 
   //Normalize the eigenvectors
   for(int i=0;i< evecs.cols(); ++i)
   {
      evecs.col(i) = evecs.col(i)/sqrt(fabs(evals(i)));
   }
 
   ff.write(outFile, evecs);
}
 
template<typename eigenArrT1, typename eigenArrT2, typename eigenArrT3>
void makeKL( eigenArrT1 & kl,
             eigenArrT2 & evecs,
             eigenArrT3 && rvecs )
{
 
   int tNims = evecs.rows();
   int tNpix = rvecs.rows();
 
   int n_modes = tNims;
 
   //Now calculate KL images
   /*
    *  KL = E^T * R  ==> C = A^T * B
    */
   math::gemm<typename eigenArrT1::Scalar>(CblasColMajor, CblasTrans, CblasTrans, n_modes, tNpix,
                              tNims, 1., evecs.data(), evecs.rows(), rvecs.data(), rvecs.rows(),
                                 0., kl.data(), kl.rows());
 
}
 
template<typename realT>
void makeFKL( const std::string & outFile,
              const std::string & coeffs,
              int N,
              int pupSize )
{
   fits::fitsFile<realT> ff;
   Eigen::Array<realT, -1, -1> evecs;
 
   ff.read(evecs, coeffs);
 
   improc::eigenCube<realT> Rims;
   sigproc::makeFourierBasis_Rect(Rims, pupSize, N, MX_FOURIER_MODIFIED);
 
   std::cout << Rims.planes() << " " << evecs.cols() << "\n";
 
   Eigen::Array<realT,-1,-1> kl;
   kl.resize( Rims.planes(), Rims.rows()*Rims.cols());
 
   std::cerr  << 1 << "\n";
   makeKL( kl, evecs, Rims.cube());
   std::cerr  << 2 << "\n";
 
   Eigen::Array<realT,-1,-1> klT = kl.transpose();
   //kl.resize(0,0);
   //Rims.resize(0,0);
   //evecs.resize(0,0);
 
   improc::eigenCube<realT> klims(klT.data(), Rims.rows(), Rims.cols(), Rims.planes());
 
   improc::eigenCube<realT> klimsR;
   klimsR.resize( klims.rows(), klims.cols(), klims.planes());
 
   for(int i=0; i< klims.planes(); ++i)
   {
      klimsR.image(i) = klims.image(klims.planes()-1-i);
   }
 
   ff.write(outFile, klimsR);
   std::cerr  << 3 << "\n";
}
 
} //namespace analysis
} //namespace AO
} //namespace mx
 
#endif //fourierCovariance_hpp