Line data    Source code

       1             : 
       2             : // -----------------------------------------------------------------------------
       3             : 
       4             : #include <cstdlib>
       5             : 
       6             : #include <math.h>
       7             : double rgauss( void );
       8             : 
       9             : #include <calanalysis/CalAnalysis/CalStats.h>
      10             : 
      11             : using namespace casacore;
      12             : using namespace casa;
      13             : 
      14             : // -----------------------------------------------------------------------------
      15             : 
      16           1 : int main( void ) {
      17             : 
      18             :   // Initialize the shape of the input cubes
      19             : 
      20           1 :   uInt uiNumPol = 2;
      21           1 :   uInt uiNumFreq = 64;
      22           1 :   uInt uiNumTime = 100;
      23             : 
      24           1 :   IPosition oShape( 3, uiNumPol, uiNumFreq, uiNumTime );
      25             : 
      26             : 
      27             :   // Initialize the input data cube (constant across frequency axis)
      28             : 
      29           1 :   Cube<Double> oData( oShape );
      30             : 
      31           3 :   for ( uInt p=0; p<uiNumPol; p++ ) {
      32         202 :     for ( uInt t=0; t<uiNumTime; t++ ) {
      33         200 :       Double dData = 10.0 * ( ((Double) rand())/((Double) RAND_MAX) - 0.5);
      34       13000 :       for ( uInt f=0; f<uiNumFreq; f++ ) {
      35       12800 :         oData.operator()(p,f,t) = dData + 0.1*rgauss();
      36             :       }
      37             :     }
      38             :   }
      39             : 
      40             : 
      41             :   // Initialize the input data error cube
      42             : 
      43           1 :   Cube<Double> oDataErr( oShape );
      44             : 
      45           3 :   for ( uInt p=0; p<uiNumPol; p++ ) {
      46         202 :     for ( uInt t=0; t<uiNumTime; t++ ) {
      47         200 :       Double dDataErr = 0.1;
      48       13000 :       for ( uInt f=0; f<uiNumFreq; f++ ) {
      49       12800 :         oDataErr.operator()(p,f,t) = dDataErr;
      50             :       }
      51             :     }
      52             :   }
      53             : 
      54             : 
      55             :   // Initialize the input flag cube
      56             : 
      57           1 :   Cube<Bool> oFlag( oShape, false );
      58             : 
      59             : 
      60             :   // Initialize the polarization (feed) abscissa
      61             : 
      62           1 :   Vector<String> oFeed( uiNumPol );
      63           1 :   oFeed[0] = "R";
      64           1 :   oFeed[1] = "L";
      65             : 
      66             : 
      67             :   // Initialize the frequency abscissa
      68             : 
      69           1 :   Vector<Double> oFreq( uiNumFreq );
      70          65 :   for ( uInt f=0; f<uiNumFreq; f++ ) oFreq[f] = f*2.0E+06 + 10.0E+09;
      71             : 
      72             : 
      73             :   // Initialize the time abscissa
      74             : 
      75           1 :   Vector<Double> oTime( uiNumTime );
      76         101 :   for ( uInt t=0; t<uiNumTime; t++ ) oTime[t] = (Double) t;
      77             : 
      78             : 
      79             :   // Initialize the user-supplied iteration axis (the polarization axis is
      80             :   // always an iteration axis, by default)
      81             : 
      82           1 :   CalStats::AXIS eAxis = CalStats::TIME;
      83             : 
      84             : 
      85             :   // Create a CalStatsReal object and get the data iterated along the
      86             :   // polarization and time axes
      87             : 
      88           1 :   CalStats oCS( oData, oDataErr, oFlag, oFeed, oFreq, oTime, eAxis );
      89             : 
      90             :   CalStats::ARG<CalStats::NONE> oArg;
      91             :   Matrix<CalStats::OUT<CalStats::NONE> >
      92           1 :       oMatrix( oCS.stats<CalStats::NONE>( oArg ) );
      93             : 
      94             : 
      95             :   // Write the data for each polarization and time axis
      96             : 
      97           3 :   for ( uInt p=0; p<uiNumPol; p++ ) {
      98         202 :     for ( uInt t=0; t<uiNumTime; t++ ) {
      99         200 :       cout << endl << flush;
     100         200 :       cout << p << ' ' << t << endl << flush;
     101         200 :       cout << oMatrix(p,t).oAxes.eAxisIterFeedID << ' '
     102         200 :            << oMatrix(p,t).oAxes.eAxisIterUserID << ' '
     103         200 :            << oMatrix(p,t).oAxes.eAxisNonIterID << ' '
     104         200 :            << endl << flush;
     105         200 :       cout << oMatrix(p,t).oAxes.sFeed << ' '
     106         200 :            << oMatrix(p,t).oAxes.dAxisIterUser << endl << flush;
     107         200 :       cout << oMatrix(p,t).oData.oAbs << endl << flush;
     108         200 :       cout << oMatrix(p,t).oData.oFlag << endl << flush;
     109             :     }
     110             :   }
     111             : 
     112             : 
     113             :   // Return 0
     114             : 
     115           1 :   return( 0 );
     116             : 
     117           1 : }
     118             : 
     119             : // -----------------------------------------------------------------------------
     120             : 
     121       12800 : double rgauss( void ) {
     122             : 
     123             :   static int iset=0;
     124             :   static double gset;
     125             :   double fac,r,v1,v2;
     126             : 
     127       12800 :   if ( iset == 0 ) {
     128             :     do {
     129        8206 :       v1 = 2.0 * (((double) rand())/((double) RAND_MAX) - 0.5);
     130        8206 :       v2 = 2.0 * (((double) rand())/((double) RAND_MAX) - 0.5);
     131        8206 :       r = v1*v1 + v2*v2;
     132        8206 :     } while ( r >= 1.0 || r == 0.0 );
     133        6400 :     fac = sqrt( -2.0*log(r)/r );
     134        6400 :     gset = v1*fac;
     135        6400 :     iset = 1;
     136        6400 :     return( v2*fac );
     137             :   } else {
     138        6400 :     iset = 0;
     139        6400 :     return( gset );
     140             :   }
     141             : 
     142             : }