[2356] | 1 | #include <iostream>
|
---|
| 2 | #include <fstream>
|
---|
| 3 |
|
---|
| 4 | #include <casa/BasicSL/String.h>
|
---|
| 5 | #include <casa/Arrays/Vector.h>
|
---|
| 6 | #include <casa/Arrays/Matrix.h>
|
---|
| 7 | #include <casa/Arrays/Cube.h>
|
---|
| 8 | #include <casa/Arrays/ArrayMath.h>
|
---|
[2361] | 9 | #include <casa/Arrays/ArrayPartMath.h>
|
---|
[2356] | 10 | #include <casa/Quanta/Quantum.h>
|
---|
| 11 | #include <casa/Quanta/QuantumHolder.h>
|
---|
| 12 | #include <casa/Utilities/CountedPtr.h>
|
---|
[2361] | 13 | #include <casa/Logging/LogIO.h>
|
---|
[2356] | 14 |
|
---|
| 15 | #include <tables/Tables/Table.h>
|
---|
[2360] | 16 | #include <tables/Tables/TableRecord.h>
|
---|
| 17 | #include <tables/Tables/ExprNode.h>
|
---|
[2356] | 18 | #include <tables/Tables/ScalarColumn.h>
|
---|
| 19 | #include <tables/Tables/ArrayColumn.h>
|
---|
| 20 |
|
---|
| 21 | #include <measures/Measures/MDirection.h>
|
---|
| 22 |
|
---|
[2364] | 23 | #include <MathUtils.h>
|
---|
| 24 |
|
---|
[2356] | 25 | #include "STGrid.h"
|
---|
| 26 |
|
---|
| 27 | using namespace std ;
|
---|
| 28 | using namespace casa ;
|
---|
| 29 | using namespace asap ;
|
---|
| 30 |
|
---|
| 31 | namespace asap {
|
---|
| 32 |
|
---|
| 33 | // constructor
|
---|
| 34 | STGrid::STGrid()
|
---|
| 35 | {
|
---|
| 36 | init() ;
|
---|
| 37 | }
|
---|
| 38 |
|
---|
| 39 | STGrid::STGrid( const string infile )
|
---|
| 40 | {
|
---|
| 41 | init() ;
|
---|
| 42 |
|
---|
| 43 | setFileIn( infile ) ;
|
---|
| 44 | }
|
---|
| 45 |
|
---|
| 46 | void STGrid::init()
|
---|
| 47 | {
|
---|
[2362] | 48 | ifno_ = -1 ;
|
---|
[2356] | 49 | nx_ = -1 ;
|
---|
| 50 | ny_ = -1 ;
|
---|
[2361] | 51 | npol_ = 0 ;
|
---|
| 52 | nchan_ = 0 ;
|
---|
| 53 | nrow_ = 0 ;
|
---|
[2356] | 54 | cellx_ = 0.0 ;
|
---|
| 55 | celly_ = 0.0 ;
|
---|
| 56 | center_ = Vector<Double> ( 2, 0.0 ) ;
|
---|
| 57 | convType_ = "BOX" ;
|
---|
[2361] | 58 | wtype_ = "UNIFORM" ;
|
---|
[2356] | 59 | convSupport_ = -1 ;
|
---|
| 60 | userSupport_ = -1 ;
|
---|
| 61 | convSampling_ = 100 ;
|
---|
| 62 | }
|
---|
| 63 |
|
---|
| 64 | void STGrid::setFileIn( const string infile )
|
---|
| 65 | {
|
---|
| 66 | String name( infile ) ;
|
---|
| 67 | if ( infile_.compare( name ) != 0 ) {
|
---|
| 68 | infile_ = String( infile ) ;
|
---|
| 69 | tab_ = Table( infile_ ) ;
|
---|
| 70 | }
|
---|
| 71 | }
|
---|
| 72 |
|
---|
[2360] | 73 | void STGrid::setPolList( vector<unsigned int> pols )
|
---|
| 74 | {
|
---|
| 75 | pollist_.assign( Vector<uInt>( pols ) ) ;
|
---|
| 76 | cout << "pollist_ = " << pollist_ << endl ;
|
---|
| 77 | }
|
---|
| 78 |
|
---|
[2364] | 79 | void STGrid::setScanList( vector<unsigned int> scans )
|
---|
| 80 | {
|
---|
| 81 | scanlist_.assign( Vector<uInt>( scans ) ) ;
|
---|
| 82 | cout << "scanlist_ = " << scanlist_ << endl ;
|
---|
| 83 | }
|
---|
| 84 |
|
---|
[2361] | 85 | void STGrid::setWeight( const string wType )
|
---|
| 86 | {
|
---|
| 87 | wtype_ = String( wType ) ;
|
---|
| 88 | wtype_.upcase() ;
|
---|
| 89 | cout << "wtype_ = " << wtype_ << endl ;
|
---|
| 90 | }
|
---|
| 91 |
|
---|
[2356] | 92 | void STGrid::defineImage( int nx,
|
---|
| 93 | int ny,
|
---|
| 94 | string scellx,
|
---|
| 95 | string scelly,
|
---|
| 96 | string scenter )
|
---|
| 97 | {
|
---|
| 98 | ROArrayColumn<Double> dirCol( tab_, "DIRECTION" ) ;
|
---|
| 99 | Matrix<Double> direction = dirCol.getColumn() ;
|
---|
| 100 | Double rmax, rmin, dmax, dmin ;
|
---|
| 101 | minMax( rmin, rmax, direction.row( 0 ) ) ;
|
---|
| 102 | minMax( dmin, dmax, direction.row( 1 ) ) ;
|
---|
| 103 |
|
---|
| 104 | Int npx = (Int)nx ;
|
---|
| 105 | Int npy = (Int)ny ;
|
---|
| 106 | String cellx( scellx ) ;
|
---|
| 107 | String celly( scelly ) ;
|
---|
| 108 | String center( scenter ) ;
|
---|
| 109 | setupGrid( npx, npy,
|
---|
| 110 | cellx, celly,
|
---|
| 111 | rmin, rmax,
|
---|
| 112 | dmin, dmax,
|
---|
| 113 | center ) ;
|
---|
| 114 | }
|
---|
| 115 |
|
---|
[2364] | 116 | void STGrid::setFunc( string convType,
|
---|
| 117 | int convSupport )
|
---|
[2356] | 118 | {
|
---|
| 119 | convType_ = String( convType ) ;
|
---|
| 120 | convType_.upcase() ;
|
---|
| 121 | userSupport_ = (Int)convSupport ;
|
---|
| 122 | }
|
---|
| 123 |
|
---|
| 124 | #define NEED_UNDERSCORES
|
---|
| 125 | #if defined(NEED_UNDERSCORES)
|
---|
| 126 | #define ggridsd ggridsd_
|
---|
| 127 | #endif
|
---|
| 128 | extern "C" {
|
---|
| 129 | void ggridsd(Double*,
|
---|
| 130 | const Complex*,
|
---|
| 131 | Int*,
|
---|
| 132 | Int*,
|
---|
| 133 | Int*,
|
---|
| 134 | const Int*,
|
---|
| 135 | const Int*,
|
---|
| 136 | const Float*,
|
---|
| 137 | Int*,
|
---|
| 138 | Int*,
|
---|
| 139 | Complex*,
|
---|
| 140 | Float*,
|
---|
| 141 | Int*,
|
---|
| 142 | Int*,
|
---|
| 143 | Int *,
|
---|
| 144 | Int *,
|
---|
| 145 | Int*,
|
---|
| 146 | Int*,
|
---|
| 147 | Float*,
|
---|
| 148 | Int*,
|
---|
| 149 | Int*,
|
---|
| 150 | Double*);
|
---|
| 151 | }
|
---|
| 152 | void STGrid::grid()
|
---|
| 153 | {
|
---|
[2361] | 154 | LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
|
---|
[2356] | 155 |
|
---|
| 156 | // retrieve data
|
---|
[2360] | 157 | Cube<Float> spectra ;
|
---|
[2356] | 158 | Matrix<Double> direction ;
|
---|
[2360] | 159 | Cube<uChar> flagtra ;
|
---|
| 160 | Matrix<uInt> rflag ;
|
---|
[2361] | 161 | Matrix<Float> weight ;
|
---|
[2364] | 162 | double t0, t1 ;
|
---|
| 163 | t0 = mathutil::gettimeofday_sec() ;
|
---|
[2362] | 164 | getData( spectra, direction, flagtra, rflag, weight ) ;
|
---|
[2364] | 165 | t1 = mathutil::gettimeofday_sec() ;
|
---|
| 166 | os << "getData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
|
---|
[2356] | 167 | IPosition sshape = spectra.shape() ;
|
---|
[2361] | 168 | //os << "spectra.shape()=" << spectra.shape() << LogIO::POST ;
|
---|
| 169 | //os << "max(spectra) = " << max(spectra) << LogIO::POST ;
|
---|
| 170 | //os << "weight = " << weight << LogIO::POST ;
|
---|
[2356] | 171 |
|
---|
| 172 | // flagtra: uChar -> Int
|
---|
| 173 | // rflag: uInt -> Int
|
---|
[2360] | 174 | Cube<Int> flagI ;
|
---|
| 175 | Matrix<Int> rflagI ;
|
---|
[2364] | 176 | t0 = mathutil::gettimeofday_sec() ;
|
---|
[2360] | 177 | toInt( &flagtra, &flagI ) ;
|
---|
| 178 | toInt( &rflag, &rflagI ) ;
|
---|
[2364] | 179 | t1 = mathutil::gettimeofday_sec() ;
|
---|
| 180 | os << "toInt: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
|
---|
[2356] | 181 |
|
---|
| 182 | // grid parameter
|
---|
[2362] | 183 | os << LogIO::DEBUGGING ;
|
---|
| 184 | os << "----------" << endl ;
|
---|
| 185 | os << "Grid parameter summary" << endl ;
|
---|
| 186 | os << " (nx,ny) = (" << nx_ << "," << ny_ << ")" << endl ;
|
---|
| 187 | os << " (cellx,celly) = (" << cellx_ << "," << celly_ << ")" << endl ;
|
---|
| 188 | os << " center = " << center_ << endl ;
|
---|
| 189 | os << "----------" << LogIO::POST ;
|
---|
| 190 | os << LogIO::NORMAL ;
|
---|
[2356] | 191 |
|
---|
[2359] | 192 | // convolution kernel
|
---|
| 193 | Vector<Float> convFunc ;
|
---|
[2364] | 194 | t0 = mathutil::gettimeofday_sec() ;
|
---|
[2359] | 195 | setConvFunc( convFunc ) ;
|
---|
[2364] | 196 | t1 = mathutil::gettimeofday_sec() ;
|
---|
| 197 | os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
|
---|
[2359] | 198 | //cout << "convSupport=" << convSupport_ << endl ;
|
---|
| 199 | //cout << "convFunc=" << convFunc << endl ;
|
---|
| 200 |
|
---|
[2356] | 201 | // world -> pixel
|
---|
| 202 | Matrix<Double> xypos( direction.shape(), 0.0 ) ;
|
---|
[2364] | 203 | t0 = mathutil::gettimeofday_sec() ;
|
---|
[2356] | 204 | toPixel( direction, xypos ) ;
|
---|
[2364] | 205 | t1 = mathutil::gettimeofday_sec() ;
|
---|
| 206 | os << "toPixel: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
|
---|
[2356] | 207 |
|
---|
[2361] | 208 | // call ggridsd
|
---|
[2356] | 209 | Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG ;
|
---|
| 210 | Double *xypos_p = xypos.getStorage( deletePos ) ;
|
---|
[2360] | 211 | Cube<Complex> dataC( spectra.shape(), 0.0 ) ;
|
---|
[2356] | 212 | setReal( dataC, spectra ) ;
|
---|
| 213 | const Complex *data_p = dataC.getStorage( deleteData ) ;
|
---|
| 214 | const Float *wgt_p = weight.getStorage( deleteWgt ) ;
|
---|
| 215 | const Int *flag_p = flagI.getStorage( deleteFlag ) ;
|
---|
| 216 | const Int *rflag_p = rflagI.getStorage( deleteFlagR ) ;
|
---|
| 217 | Float *conv_p = convFunc.getStorage( deleteConv ) ;
|
---|
[2359] | 218 | // Extend grid plane with convSupport_
|
---|
[2364] | 219 | Int gnx = nx_ ;
|
---|
| 220 | Int gny = ny_ ;
|
---|
| 221 | // Int gnx = nx_+convSupport_*2 ;
|
---|
| 222 | // Int gny = ny_+convSupport_*2 ;
|
---|
[2361] | 223 | IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
|
---|
[2356] | 224 | Array<Complex> gdataArrC( gshape, 0.0 ) ;
|
---|
| 225 | Array<Float> gwgtArr( gshape, 0.0 ) ;
|
---|
| 226 | Complex *gdata_p = gdataArrC.getStorage( deleteDataG ) ;
|
---|
| 227 | Float *wdata_p = gwgtArr.getStorage( deleteWgtG ) ;
|
---|
| 228 | Int idopsf = 0 ;
|
---|
[2361] | 229 | Int *chanMap = new Int[nchan_] ;
|
---|
[2356] | 230 | {
|
---|
| 231 | Int *work_p = chanMap ;
|
---|
[2361] | 232 | for ( Int i = 0 ; i < nchan_ ; i++ ) {
|
---|
[2356] | 233 | *work_p = i ;
|
---|
| 234 | work_p++ ;
|
---|
| 235 | }
|
---|
| 236 | }
|
---|
[2361] | 237 | Int *polMap = new Int[npol_] ;
|
---|
[2356] | 238 | {
|
---|
| 239 | Int *work_p = polMap ;
|
---|
[2361] | 240 | for ( Int i = 0 ; i < npol_ ; i++ ) {
|
---|
[2356] | 241 | *work_p = i ;
|
---|
| 242 | work_p++ ;
|
---|
| 243 | }
|
---|
| 244 | }
|
---|
[2361] | 245 | Double *sumw_p = new Double[npol_*nchan_] ;
|
---|
[2356] | 246 | {
|
---|
| 247 | Double *work_p = sumw_p ;
|
---|
[2361] | 248 | for ( Int i = 0 ; i < npol_*nchan_ ; i++ ) {
|
---|
[2356] | 249 | *work_p = 0.0 ;
|
---|
| 250 | work_p++ ;
|
---|
| 251 | }
|
---|
| 252 | }
|
---|
[2364] | 253 | t0 = mathutil::gettimeofday_sec() ;
|
---|
| 254 | Int irow = -1 ;
|
---|
[2356] | 255 | ggridsd( xypos_p,
|
---|
| 256 | data_p,
|
---|
[2361] | 257 | &npol_,
|
---|
| 258 | &nchan_,
|
---|
[2356] | 259 | &idopsf,
|
---|
| 260 | flag_p,
|
---|
| 261 | rflag_p,
|
---|
| 262 | wgt_p,
|
---|
[2361] | 263 | &nrow_,
|
---|
[2356] | 264 | &irow,
|
---|
| 265 | gdata_p,
|
---|
| 266 | wdata_p,
|
---|
[2359] | 267 | &gnx,
|
---|
| 268 | &gny,
|
---|
[2361] | 269 | &npol_,
|
---|
| 270 | &nchan_,
|
---|
[2356] | 271 | &convSupport_,
|
---|
| 272 | &convSampling_,
|
---|
| 273 | conv_p,
|
---|
| 274 | chanMap,
|
---|
| 275 | polMap,
|
---|
| 276 | sumw_p ) ;
|
---|
[2364] | 277 | t1 = mathutil::gettimeofday_sec() ;
|
---|
| 278 | os << "ggridsd: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
|
---|
[2356] | 279 | xypos.putStorage( xypos_p, deletePos ) ;
|
---|
| 280 | dataC.freeStorage( data_p, deleteData ) ;
|
---|
| 281 | weight.freeStorage( wgt_p, deleteWgt ) ;
|
---|
| 282 | flagI.freeStorage( flag_p, deleteFlag ) ;
|
---|
| 283 | rflagI.freeStorage( rflag_p, deleteFlagR ) ;
|
---|
| 284 | convFunc.putStorage( conv_p, deleteConv ) ;
|
---|
| 285 | delete polMap ;
|
---|
| 286 | delete chanMap ;
|
---|
| 287 | gdataArrC.putStorage( gdata_p, deleteDataG ) ;
|
---|
| 288 | gwgtArr.putStorage( wdata_p, deleteWgtG ) ;
|
---|
| 289 | Array<Float> gdataArr = real( gdataArrC ) ;
|
---|
[2368] | 290 | setData( data_, gdataArr, gwgtArr ) ;
|
---|
[2361] | 291 | //Matrix<Double> sumWeight( IPosition( 2, npol_, nchan_ ), sumw_p, TAKE_OVER ) ;
|
---|
| 292 | delete sumw_p ;
|
---|
[2356] | 293 | //cout << "sumWeight = " << sumWeight << endl ;
|
---|
[2364] | 294 | // os << "gdataArr = " << gdataArr << LogIO::POST ;
|
---|
| 295 | // os << "gwgtArr = " << gwgtArr << LogIO::POST ;
|
---|
| 296 | // os << "data_ " << data_ << LogIO::POST ;
|
---|
[2356] | 297 | }
|
---|
| 298 |
|
---|
[2368] | 299 | void STGrid::setData( Array<Float> &data,
|
---|
| 300 | Array<Float> &gdata,
|
---|
| 301 | Array<Float> &gwgt )
|
---|
| 302 | {
|
---|
| 303 | LogIO os( LogOrigin("STGrid","setData",WHERE) ) ;
|
---|
| 304 | double t0, t1 ;
|
---|
| 305 | t0 = mathutil::gettimeofday_sec() ;
|
---|
[2371] | 306 | data.resize( gdata.shape() ) ;
|
---|
[2368] | 307 | uInt len = data.nelements() ;
|
---|
| 308 | Float *w0_p ;
|
---|
| 309 | const Float *w1_p, *w2_p ;
|
---|
| 310 | Bool b0, b1, b2 ;
|
---|
| 311 | Float *data_p = data.getStorage( b0 ) ;
|
---|
| 312 | const Float *gdata_p = gdata.getStorage( b1 ) ;
|
---|
| 313 | const Float *gwgt_p = gwgt.getStorage( b2 ) ;
|
---|
| 314 | w0_p = data_p ;
|
---|
| 315 | w1_p = gdata_p ;
|
---|
| 316 | w2_p = gwgt_p ;
|
---|
| 317 | for ( uInt i = 0 ; i < len ; i++ ) {
|
---|
| 318 | *w0_p = (*w2_p > 0.0) ? (*w1_p / *w2_p) : 0.0 ;
|
---|
| 319 | w0_p++ ;
|
---|
| 320 | w1_p++ ;
|
---|
| 321 | w2_p++ ;
|
---|
| 322 | }
|
---|
| 323 | data.putStorage( data_p, b0 ) ;
|
---|
| 324 | gdata.freeStorage( gdata_p, b1 ) ;
|
---|
| 325 | gwgt.freeStorage( gwgt_p, b2 ) ;
|
---|
| 326 | t1 = mathutil::gettimeofday_sec() ;
|
---|
| 327 | os << "setData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
|
---|
| 328 | }
|
---|
| 329 |
|
---|
[2356] | 330 | void STGrid::setupGrid( Int &nx,
|
---|
| 331 | Int &ny,
|
---|
| 332 | String &cellx,
|
---|
| 333 | String &celly,
|
---|
| 334 | Double &xmin,
|
---|
| 335 | Double &xmax,
|
---|
| 336 | Double &ymin,
|
---|
| 337 | Double &ymax,
|
---|
| 338 | String ¢er )
|
---|
| 339 | {
|
---|
| 340 | //cout << "nx=" << nx << ", ny=" << ny << endl ;
|
---|
[2359] | 341 |
|
---|
| 342 | // center position
|
---|
| 343 | if ( center.size() == 0 ) {
|
---|
| 344 | center_(0) = 0.5 * ( xmin + xmax ) ;
|
---|
| 345 | center_(1) = 0.5 * ( ymin + ymax ) ;
|
---|
| 346 | }
|
---|
| 347 | else {
|
---|
| 348 | String::size_type pos0 = center.find( " " ) ;
|
---|
| 349 | if ( pos0 == String::npos ) {
|
---|
| 350 | throw AipsError( "bad string format in parameter center" ) ;
|
---|
| 351 | }
|
---|
| 352 | String::size_type pos1 = center.find( " ", pos0+1 ) ;
|
---|
| 353 | String typestr, xstr, ystr ;
|
---|
| 354 | if ( pos1 != String::npos ) {
|
---|
| 355 | typestr = center.substr( 0, pos0 ) ;
|
---|
| 356 | xstr = center.substr( pos0+1, pos1-pos0 ) ;
|
---|
| 357 | ystr = center.substr( pos1+1 ) ;
|
---|
| 358 | // todo: convert to J2000 (or direction ref for DIRECTION column)
|
---|
| 359 | }
|
---|
| 360 | else {
|
---|
| 361 | typestr = "J2000" ;
|
---|
| 362 | xstr = center.substr( 0, pos0 ) ;
|
---|
| 363 | ystr = center.substr( pos0+1 ) ;
|
---|
| 364 | }
|
---|
| 365 | QuantumHolder qh ;
|
---|
| 366 | String err ;
|
---|
| 367 | qh.fromString( err, xstr ) ;
|
---|
| 368 | Quantum<Double> xcen = qh.asQuantumDouble() ;
|
---|
| 369 | qh.fromString( err, ystr ) ;
|
---|
| 370 | Quantum<Double> ycen = qh.asQuantumDouble() ;
|
---|
| 371 | center_(0) = xcen.getValue( "rad" ) ;
|
---|
| 372 | center_(1) = ycen.getValue( "rad" ) ;
|
---|
| 373 | }
|
---|
| 374 |
|
---|
| 375 |
|
---|
| 376 | //Double wx = xmax - xmin ;
|
---|
| 377 | //Double wy = ymax - ymin ;
|
---|
| 378 | Double wx = max( abs(xmax-center_(0)), abs(xmin-center_(0)) ) * 2 ;
|
---|
| 379 | Double wy = max( abs(ymax-center_(1)), abs(ymin-center_(1)) ) * 2 ;
|
---|
| 380 | // take 10% margin
|
---|
[2356] | 381 | wx *= 1.10 ;
|
---|
| 382 | wy *= 1.10 ;
|
---|
| 383 | Quantum<Double> qcellx ;
|
---|
| 384 | Quantum<Double> qcelly ;
|
---|
| 385 | nx_ = nx ;
|
---|
| 386 | ny_ = ny ;
|
---|
| 387 | if ( nx < 0 && ny > 0 ) {
|
---|
| 388 | nx_ = ny ;
|
---|
| 389 | ny_ = ny ;
|
---|
| 390 | }
|
---|
| 391 | if ( ny < 0 && nx > 0 ) {
|
---|
| 392 | nx_ = nx ;
|
---|
| 393 | ny_ = nx ;
|
---|
| 394 | }
|
---|
| 395 | //cout << "nx_ = " << nx_ << ", ny_ = " << ny_ << endl ;
|
---|
| 396 | if ( cellx.size() != 0 && celly.size() != 0 ) {
|
---|
| 397 | readQuantity( qcellx, cellx ) ;
|
---|
| 398 | readQuantity( qcelly, celly ) ;
|
---|
| 399 | }
|
---|
| 400 | else if ( celly.size() != 0 ) {
|
---|
| 401 | cout << "Using celly to x-axis..." << endl ;
|
---|
| 402 | readQuantity( qcelly, celly ) ;
|
---|
| 403 | qcellx = qcelly ;
|
---|
| 404 | }
|
---|
| 405 | else if ( cellx.size() != 0 ) {
|
---|
| 406 | cout << "Using cellx to y-axis..." << endl ;
|
---|
| 407 | readQuantity( qcellx, cellx ) ;
|
---|
| 408 | qcelly = qcellx ;
|
---|
| 409 | }
|
---|
| 410 | else {
|
---|
| 411 | if ( nx_ < 0 ) {
|
---|
| 412 | cout << "No user preference in grid setting. Using default..." << endl ;
|
---|
| 413 | readQuantity( qcellx, "1.0arcmin" ) ;
|
---|
| 414 | qcelly = qcellx ;
|
---|
| 415 | }
|
---|
| 416 | else {
|
---|
| 417 | qcellx = Quantum<Double>( wx/nx_, "rad" ) ;
|
---|
| 418 | qcelly = Quantum<Double>( wy/ny_, "rad" ) ;
|
---|
| 419 | }
|
---|
| 420 | }
|
---|
| 421 | cellx_ = qcellx.getValue( "rad" ) ;
|
---|
| 422 | celly_ = qcelly.getValue( "rad" ) ;
|
---|
| 423 | if ( nx_ < 0 ) {
|
---|
| 424 | nx_ = Int( ceil( wx/cellx_ ) ) ;
|
---|
| 425 | ny_ = Int( ceil( wy/celly_ ) ) ;
|
---|
| 426 | }
|
---|
| 427 | }
|
---|
| 428 |
|
---|
[2362] | 429 | void STGrid::selectData( Table &tab )
|
---|
| 430 | {
|
---|
| 431 | Int ifno = ifno_ ;
|
---|
| 432 | Table taborg( infile_ ) ;
|
---|
| 433 | if ( ifno == -1 ) {
|
---|
[2368] | 434 | LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;
|
---|
[2362] | 435 | // os << LogIO::SEVERE
|
---|
| 436 | // << "Please set IFNO before actual gridding"
|
---|
| 437 | // << LogIO::EXCEPTION ;
|
---|
| 438 | ROScalarColumn<uInt> ifnoCol( taborg, "IFNO" ) ;
|
---|
| 439 | ifno = ifnoCol( 0 ) ;
|
---|
| 440 | os << LogIO::WARN
|
---|
| 441 | << "IFNO is not given. Using default IFNO: " << ifno << LogIO::POST ;
|
---|
| 442 | }
|
---|
[2364] | 443 | // tab = taborg( taborg.col("IFNO") == ifno ) ;
|
---|
| 444 | TableExprNode node ;
|
---|
| 445 | node = taborg.col("IFNO") == ifno ;
|
---|
| 446 | if ( scanlist_.size() > 0 ) {
|
---|
| 447 | node = node && taborg.col("SCANNO").in( scanlist_ ) ;
|
---|
| 448 | }
|
---|
| 449 | tab = taborg( node ) ;
|
---|
[2362] | 450 | if ( tab.nrow() == 0 ) {
|
---|
[2368] | 451 | LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;
|
---|
[2362] | 452 | os << LogIO::SEVERE
|
---|
[2365] | 453 | << "No corresponding rows for given selection: IFNO " << ifno
|
---|
| 454 | << " SCANNO " << scanlist_
|
---|
[2362] | 455 | << LogIO::EXCEPTION ;
|
---|
| 456 | }
|
---|
| 457 | }
|
---|
| 458 |
|
---|
| 459 | void STGrid::getData( Cube<Float> &spectra,
|
---|
[2356] | 460 | Matrix<Double> &direction,
|
---|
[2360] | 461 | Cube<uChar> &flagtra,
|
---|
[2361] | 462 | Matrix<uInt> &rflag,
|
---|
| 463 | Matrix<Float> &weight )
|
---|
[2356] | 464 | {
|
---|
[2362] | 465 | Table tab ;
|
---|
| 466 | selectData( tab ) ;
|
---|
[2371] | 467 | updatePolList( tab ) ;
|
---|
[2361] | 468 | // cout << "npol_ = " << npol_ << endl ;
|
---|
| 469 | // cout << "nchan_ = " << nchan_ << endl ;
|
---|
| 470 | // cout << "nrow_ = " << nrow_ << endl ;
|
---|
| 471 | spectra.resize( npol_, nchan_, nrow_ ) ;
|
---|
| 472 | flagtra.resize( npol_, nchan_, nrow_ ) ;
|
---|
| 473 | rflag.resize( npol_, nrow_ ) ;
|
---|
| 474 | Cube<Float> tsys( npol_, nchan_, nrow_ ) ;
|
---|
| 475 | Matrix<Double> tint( npol_, nrow_ ) ;
|
---|
[2371] | 476 | // boolean for pointer access
|
---|
| 477 | Bool bsp, bfl, bfr, bts, bti ;
|
---|
| 478 | // pointer to the data
|
---|
| 479 | Float *sp_p = spectra.getStorage( bsp ) ;
|
---|
| 480 | uChar *fl_p = flagtra.getStorage( bfl ) ;
|
---|
| 481 | uInt *fr_p = rflag.getStorage( bfr ) ;
|
---|
| 482 | Float *ts_p = tsys.getStorage( bts ) ;
|
---|
| 483 | Double *ti_p = tint.getStorage( bti ) ;
|
---|
| 484 | // working pointer
|
---|
| 485 | Float *wsp_p = sp_p ;
|
---|
| 486 | uChar *wfl_p = fl_p ;
|
---|
| 487 | uInt *wfr_p = fr_p ;
|
---|
| 488 | Float *wts_p = ts_p ;
|
---|
| 489 | Double *wti_p = ti_p ;
|
---|
| 490 | uInt len = nchan_ * nrow_ ;
|
---|
| 491 | IPosition mshape( 2, nchan_, nrow_ ) ;
|
---|
| 492 | IPosition vshape( 1, nrow_ ) ;
|
---|
[2361] | 493 | for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
|
---|
[2360] | 494 | Table subt = tab( tab.col("POLNO") == pollist_[ipol] ) ;
|
---|
| 495 | ROArrayColumn<Float> spectraCol( subt, "SPECTRA" ) ;
|
---|
| 496 | ROArrayColumn<Double> directionCol( subt, "DIRECTION" ) ;
|
---|
| 497 | ROArrayColumn<uChar> flagtraCol( subt, "FLAGTRA" ) ;
|
---|
| 498 | ROScalarColumn<uInt> rflagCol( subt, "FLAGROW" ) ;
|
---|
[2361] | 499 | ROArrayColumn<Float> tsysCol( subt, "TSYS" ) ;
|
---|
| 500 | ROScalarColumn<Double> tintCol( subt, "INTERVAL" ) ;
|
---|
[2371] | 501 | Matrix<Float> spSlice( mshape, wsp_p, SHARE ) ;
|
---|
| 502 | Matrix<uChar> flSlice( mshape, wfl_p, SHARE ) ;
|
---|
| 503 | Vector<uInt> frSlice( vshape, wfr_p, SHARE ) ;
|
---|
| 504 | spectraCol.getColumn( spSlice ) ;
|
---|
| 505 | flagtraCol.getColumn( flSlice ) ;
|
---|
| 506 | rflagCol.getColumn( frSlice ) ;
|
---|
[2360] | 507 | if ( ipol == 0 )
|
---|
| 508 | directionCol.getColumn( direction ) ;
|
---|
[2371] | 509 | Vector<Float> tmpF = tsysCol( 0 ) ;
|
---|
| 510 | Vector<Double> tmpD( vshape, wti_p, SHARE ) ;
|
---|
| 511 | Matrix<Float> tsSlice( mshape, wts_p, SHARE ) ;
|
---|
| 512 | if ( tmpF.nelements() == (uInt)nchan_ ) {
|
---|
| 513 | tsysCol.getColumn( tsSlice ) ;
|
---|
[2361] | 514 | }
|
---|
| 515 | else {
|
---|
[2371] | 516 | tsSlice = tmpF( 0 ) ;
|
---|
[2361] | 517 | }
|
---|
| 518 | tintCol.getColumn( tmpD ) ;
|
---|
[2371] | 519 |
|
---|
| 520 | wsp_p += len ;
|
---|
| 521 | wfl_p += len ;
|
---|
| 522 | wfr_p += nrow_ ;
|
---|
| 523 | wts_p += len ;
|
---|
| 524 | wti_p += nrow_ ;
|
---|
[2360] | 525 | }
|
---|
[2371] | 526 | spectra.putStorage( sp_p, bsp ) ;
|
---|
| 527 | flagtra.putStorage( fl_p, bfl ) ;
|
---|
| 528 | rflag.putStorage( fr_p, bfr ) ;
|
---|
| 529 | tsys.putStorage( ts_p, bts ) ;
|
---|
| 530 | tint.putStorage( ti_p, bti ) ;
|
---|
[2361] | 531 |
|
---|
| 532 | getWeight( weight, tsys, tint ) ;
|
---|
[2356] | 533 | }
|
---|
| 534 |
|
---|
[2371] | 535 | void STGrid::updatePolList( Table &tab )
|
---|
| 536 | {
|
---|
| 537 | ROScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
|
---|
| 538 | Vector<uInt> pols = polnoCol.getColumn() ;
|
---|
| 539 | Vector<uInt> pollistOrg ;
|
---|
| 540 | uInt npolOrg = 0 ;
|
---|
| 541 | uInt polno ;
|
---|
| 542 | for ( uInt i = 0 ; i < polnoCol.nrow() ; i++ ) {
|
---|
| 543 | //polno = polnoCol( i ) ;
|
---|
| 544 | polno = pols( i ) ;
|
---|
| 545 | if ( allNE( pollistOrg, polno ) ) {
|
---|
| 546 | pollistOrg.resize( npolOrg+1, True ) ;
|
---|
| 547 | pollistOrg[npolOrg] = polno ;
|
---|
| 548 | npolOrg++ ;
|
---|
| 549 | }
|
---|
| 550 | }
|
---|
| 551 | if ( pollist_.size() == 0 )
|
---|
| 552 | pollist_ = pollistOrg ;
|
---|
| 553 | else {
|
---|
| 554 | Vector<uInt> newlist ;
|
---|
| 555 | uInt newsize = 0 ;
|
---|
| 556 | for ( uInt i = 0 ; i < pollist_.size() ; i++ ) {
|
---|
| 557 | if ( anyEQ( pollistOrg, pollist_[i] ) ) {
|
---|
| 558 | newlist.resize( newsize+1, True ) ;
|
---|
| 559 | newlist[newsize] = pollist_[i] ;
|
---|
| 560 | newsize++ ;
|
---|
| 561 | }
|
---|
| 562 | }
|
---|
| 563 | pollist_.assign( newlist ) ;
|
---|
| 564 | }
|
---|
| 565 | npol_ = pollist_.size() ;
|
---|
| 566 | if ( npol_ == 0 ) {
|
---|
| 567 | LogIO os( LogOrigin("STGrid","updatePolList",WHERE) ) ;
|
---|
| 568 | os << LogIO::SEVERE << "Empty pollist" << LogIO::EXCEPTION ;
|
---|
| 569 | }
|
---|
| 570 | nrow_ = tab.nrow() / npolOrg ;
|
---|
| 571 | ROArrayColumn<uChar> tmpCol( tab, "FLAGTRA" ) ;
|
---|
| 572 | nchan_ = tmpCol( 0 ).nelements() ;
|
---|
| 573 | // LogIO os( LogOrigin("STGrid","updatePolList",WHERE) ) ;
|
---|
| 574 | // os << "npol_ = " << npol_ << "(" << pollist_ << ")" << endl
|
---|
| 575 | // << "nchan_ = " << nchan_ << endl
|
---|
| 576 | // << "nrow_ = " << nrow_ << LogIO::POST ;
|
---|
| 577 | }
|
---|
| 578 |
|
---|
[2361] | 579 | void STGrid::getWeight( Matrix<Float> &w,
|
---|
| 580 | Cube<Float> &tsys,
|
---|
| 581 | Matrix<Double> &tint )
|
---|
| 582 | {
|
---|
[2368] | 583 | LogIO os( LogOrigin("STGrid","getWeight",WHERE) ) ;
|
---|
| 584 | double t0, t1 ;
|
---|
| 585 | t0 = mathutil::gettimeofday_sec() ;
|
---|
[2361] | 586 | // resize
|
---|
| 587 | w.resize( nchan_, nrow_ ) ;
|
---|
| 588 |
|
---|
| 589 | // set weight
|
---|
| 590 | Bool warn = False ;
|
---|
| 591 | if ( wtype_.compare( "UNIFORM" ) == 0 ) {
|
---|
[2368] | 592 | w = 1.0 ;
|
---|
[2361] | 593 | }
|
---|
| 594 | else if ( wtype_.compare( "TINT" ) == 0 ) {
|
---|
| 595 | if ( npol_ > 1 ) warn = True ;
|
---|
[2369] | 596 | Bool b0, b1 ;
|
---|
| 597 | Float *w_p = w.getStorage( b0 ) ;
|
---|
| 598 | Float *w0_p = w_p ;
|
---|
| 599 | const Double *ti_p = tint.getStorage( b1 ) ;
|
---|
| 600 | const Double *w1_p = ti_p ;
|
---|
[2361] | 601 | for ( Int irow = 0 ; irow < nrow_ ; irow++ ) {
|
---|
[2369] | 602 | Float val = (Float)(polMean( w1_p )) ;
|
---|
| 603 | for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
|
---|
| 604 | *w0_p = val ;
|
---|
| 605 | w0_p++ ;
|
---|
| 606 | }
|
---|
[2361] | 607 | }
|
---|
[2369] | 608 | w.putStorage( w_p, b0 ) ;
|
---|
| 609 | tint.freeStorage( ti_p, b1 ) ;
|
---|
[2361] | 610 | }
|
---|
| 611 | else if ( wtype_.compare( "TSYS" ) == 0 ) {
|
---|
| 612 | if ( npol_ > 1 ) warn = True ;
|
---|
[2369] | 613 | Bool b0, b1 ;
|
---|
| 614 | Float *w_p = w.getStorage( b0 ) ;
|
---|
| 615 | Float *w0_p = w_p ;
|
---|
| 616 | const Float *ts_p = tsys.getStorage( b1 ) ;
|
---|
| 617 | const Float *w1_p = ts_p ;
|
---|
[2361] | 618 | for ( Int irow = 0 ; irow < nrow_ ; irow++ ) {
|
---|
| 619 | for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
|
---|
[2369] | 620 | Float val = polMean( w1_p ) ;
|
---|
| 621 | *w0_p = 1.0 / ( val * val ) ;
|
---|
| 622 | w0_p++ ;
|
---|
[2361] | 623 | }
|
---|
| 624 | }
|
---|
[2369] | 625 | w.putStorage( w_p, b0 ) ;
|
---|
| 626 | tsys.freeStorage( ts_p, b1 ) ;
|
---|
[2361] | 627 | }
|
---|
| 628 | else if ( wtype_.compare( "TINTSYS" ) == 0 ) {
|
---|
| 629 | if ( npol_ > 1 ) warn = True ;
|
---|
[2369] | 630 | Bool b0, b1, b2 ;
|
---|
| 631 | Float *w_p = w.getStorage( b0 ) ;
|
---|
| 632 | Float *w0_p = w_p ;
|
---|
| 633 | const Double *ti_p = tint.getStorage( b1 ) ;
|
---|
| 634 | const Double *w1_p = ti_p ;
|
---|
| 635 | const Float *ts_p = tsys.getStorage( b2 ) ;
|
---|
| 636 | const Float *w2_p = ts_p ;
|
---|
[2361] | 637 | for ( Int irow = 0 ; irow < nrow_ ; irow++ ) {
|
---|
[2369] | 638 | Float interval = (Float)(polMean( w1_p )) ;
|
---|
[2361] | 639 | for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
|
---|
[2369] | 640 | Float temp = polMean( w2_p ) ;
|
---|
| 641 | *w0_p = interval / ( temp * temp ) ;
|
---|
| 642 | w0_p++ ;
|
---|
[2361] | 643 | }
|
---|
| 644 | }
|
---|
[2369] | 645 | w.putStorage( w_p, b0 ) ;
|
---|
| 646 | tint.freeStorage( ti_p, b1 ) ;
|
---|
| 647 | tsys.freeStorage( ts_p, b2 ) ;
|
---|
[2361] | 648 | }
|
---|
| 649 | else {
|
---|
[2368] | 650 | //LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
|
---|
| 651 | os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
|
---|
| 652 | w = 1.0 ;
|
---|
[2361] | 653 | }
|
---|
| 654 |
|
---|
| 655 | if ( npol_ > 1 ) {
|
---|
[2368] | 656 | //LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
|
---|
[2362] | 657 | os << LogIO::WARN << "STGrid doesn't support assigning polarization-dependent weight. Use averaged weight over polarization." << LogIO::POST ;
|
---|
[2361] | 658 | }
|
---|
[2368] | 659 | t1 = mathutil::gettimeofday_sec() ;
|
---|
| 660 | os << "getWeight: elapsed time is " << t1-t0 << " sec" << LogIO::POST ;
|
---|
[2361] | 661 | }
|
---|
| 662 |
|
---|
[2369] | 663 | Float STGrid::polMean( const Float *p )
|
---|
| 664 | {
|
---|
| 665 | Float v = 0.0 ;
|
---|
| 666 | for ( Int i = 0 ; i < npol_ ; i++ ) {
|
---|
| 667 | v += *p ;
|
---|
| 668 | p++ ;
|
---|
| 669 | }
|
---|
| 670 | v /= npol_ ;
|
---|
| 671 | return v ;
|
---|
| 672 | }
|
---|
| 673 |
|
---|
| 674 | Double STGrid::polMean( const Double *p )
|
---|
| 675 | {
|
---|
| 676 | Double v = 0.0 ;
|
---|
| 677 | for ( Int i = 0 ; i < npol_ ; i++ ) {
|
---|
| 678 | v += *p ;
|
---|
| 679 | p++ ;
|
---|
| 680 | }
|
---|
| 681 | v /= npol_ ;
|
---|
| 682 | return v ;
|
---|
| 683 | }
|
---|
| 684 |
|
---|
[2360] | 685 | void STGrid::toInt( Array<uChar> *u, Array<Int> *v )
|
---|
[2356] | 686 | {
|
---|
[2360] | 687 | uInt len = u->nelements() ;
|
---|
[2356] | 688 | Int *int_p = new Int[len] ;
|
---|
| 689 | Bool deleteIt ;
|
---|
[2360] | 690 | const uChar *data_p = u->getStorage( deleteIt ) ;
|
---|
[2356] | 691 | Int *i_p = int_p ;
|
---|
| 692 | const uChar *u_p = data_p ;
|
---|
| 693 | for ( uInt i = 0 ; i < len ; i++ ) {
|
---|
| 694 | *i_p = ( *u_p == 0 ) ? 0 : 1 ;
|
---|
| 695 | i_p++ ;
|
---|
| 696 | u_p++ ;
|
---|
| 697 | }
|
---|
[2360] | 698 | u->freeStorage( data_p, deleteIt ) ;
|
---|
| 699 | v->takeStorage( u->shape(), int_p, TAKE_OVER ) ;
|
---|
[2356] | 700 | }
|
---|
| 701 |
|
---|
[2360] | 702 | void STGrid::toInt( Array<uInt> *u, Array<Int> *v )
|
---|
[2356] | 703 | {
|
---|
[2360] | 704 | uInt len = u->nelements() ;
|
---|
[2356] | 705 | Int *int_p = new Int[len] ;
|
---|
| 706 | Bool deleteIt ;
|
---|
[2360] | 707 | const uInt *data_p = u->getStorage( deleteIt ) ;
|
---|
[2356] | 708 | Int *i_p = int_p ;
|
---|
| 709 | const uInt *u_p = data_p ;
|
---|
| 710 | for ( uInt i = 0 ; i < len ; i++ ) {
|
---|
| 711 | *i_p = ( *u_p == 0 ) ? 0 : 1 ;
|
---|
| 712 | i_p++ ;
|
---|
| 713 | u_p++ ;
|
---|
| 714 | }
|
---|
[2360] | 715 | u->freeStorage( data_p, deleteIt ) ;
|
---|
| 716 | v->takeStorage( u->shape(), int_p, TAKE_OVER ) ;
|
---|
[2356] | 717 | }
|
---|
| 718 |
|
---|
| 719 | void STGrid::toPixel( Matrix<Double> &world, Matrix<Double> &pixel )
|
---|
| 720 | {
|
---|
[2359] | 721 | // gridding will be done on (nx_+2*convSupport_) x (ny_+2*convSupport_)
|
---|
| 722 | // grid plane to avoid unexpected behavior on grid edge
|
---|
[2356] | 723 | Vector<Double> pixc( 2 ) ;
|
---|
[2364] | 724 | pixc(0) = Double( nx_-1 ) * 0.5 ;
|
---|
| 725 | pixc(1) = Double( ny_-1 ) * 0.5 ;
|
---|
| 726 | // pixc(0) = Double( nx_+2*convSupport_-1 ) * 0.5 ;
|
---|
| 727 | // pixc(1) = Double( ny_+2*convSupport_-1 ) * 0.5 ;
|
---|
[2356] | 728 | uInt nrow = world.shape()[1] ;
|
---|
| 729 | Vector<Double> cell( 2 ) ;
|
---|
| 730 | cell(0) = cellx_ ;
|
---|
| 731 | cell(1) = celly_ ;
|
---|
| 732 | for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
|
---|
| 733 | for ( uInt i = 0 ; i < 2 ; i++ ) {
|
---|
| 734 | pixel( i, irow ) = pixc(i) + ( world(i, irow) - center_(i) ) / cell(i) ;
|
---|
| 735 | }
|
---|
| 736 | }
|
---|
[2364] | 737 | // String gridfile = "grid."+convType_+"."+String::toString(convSupport_)+".dat" ;
|
---|
| 738 | // ofstream ofs( gridfile.c_str(), ios::out ) ;
|
---|
| 739 | // ofs << "center " << center_(0) << " " << pixc(0)
|
---|
| 740 | // << " " << center_(1) << " " << pixc(1) << endl ;
|
---|
| 741 | // for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
|
---|
| 742 | // ofs << irow ;
|
---|
| 743 | // for ( uInt i = 0 ; i < 2 ; i++ ) {
|
---|
| 744 | // ofs << " " << world(i, irow) << " " << pixel(i, irow) ;
|
---|
| 745 | // }
|
---|
| 746 | // ofs << endl ;
|
---|
| 747 | // }
|
---|
| 748 | // ofs.close() ;
|
---|
[2356] | 749 | }
|
---|
| 750 |
|
---|
| 751 | void STGrid::boxFunc( Vector<Float> &convFunc, Int &convSize )
|
---|
| 752 | {
|
---|
| 753 | convFunc = 0.0 ;
|
---|
| 754 | for ( Int i = 0 ; i < convSize/2 ; i++ )
|
---|
| 755 | convFunc(i) = 1.0 ;
|
---|
| 756 | }
|
---|
| 757 |
|
---|
| 758 | #define NEED_UNDERSCORES
|
---|
| 759 | #if defined(NEED_UNDERSCORES)
|
---|
| 760 | #define grdsf grdsf_
|
---|
| 761 | #endif
|
---|
| 762 | extern "C" {
|
---|
| 763 | void grdsf(Double*, Double*);
|
---|
| 764 | }
|
---|
| 765 | void STGrid::spheroidalFunc( Vector<Float> &convFunc )
|
---|
| 766 | {
|
---|
| 767 | convFunc = 0.0 ;
|
---|
| 768 | for ( Int i = 0 ; i < convSampling_*convSupport_ ; i++ ) {
|
---|
| 769 | Double nu = Double(i) / Double(convSupport_*convSampling_) ;
|
---|
| 770 | Double val ;
|
---|
| 771 | grdsf( &nu, &val ) ;
|
---|
| 772 | convFunc(i) = ( 1.0 - nu * nu ) * val ;
|
---|
| 773 | }
|
---|
| 774 | }
|
---|
| 775 |
|
---|
| 776 | void STGrid::gaussFunc( Vector<Float> &convFunc )
|
---|
| 777 | {
|
---|
| 778 | convFunc = 0.0 ;
|
---|
[2363] | 779 | // HWHM of the Gaussian is convSupport_ / 4
|
---|
| 780 | // To take into account Gaussian tail, kernel cutoff is set to 4 * HWHM
|
---|
| 781 | Int len = convSampling_ * convSupport_ ;
|
---|
| 782 | Double hwhm = len * 0.25 ;
|
---|
| 783 | for ( Int i = 0 ; i < len ; i++ ) {
|
---|
[2356] | 784 | Double val = Double(i) / hwhm ;
|
---|
| 785 | convFunc(i) = exp( -log(2)*val*val ) ;
|
---|
| 786 | }
|
---|
| 787 | }
|
---|
| 788 |
|
---|
| 789 | void STGrid::pbFunc( Vector<Float> &convFunc )
|
---|
| 790 | {
|
---|
| 791 | convFunc = 0.0 ;
|
---|
| 792 | }
|
---|
| 793 |
|
---|
| 794 | void STGrid::setConvFunc( Vector<Float> &convFunc )
|
---|
| 795 | {
|
---|
| 796 | convSupport_ = userSupport_ ;
|
---|
| 797 | if ( convType_ == "BOX" ) {
|
---|
| 798 | if ( convSupport_ < 0 )
|
---|
| 799 | convSupport_ = 0 ;
|
---|
| 800 | Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
|
---|
| 801 | convFunc.resize( convSize ) ;
|
---|
| 802 | boxFunc( convFunc, convSize ) ;
|
---|
| 803 | }
|
---|
| 804 | else if ( convType_ == "SF" ) {
|
---|
| 805 | if ( convSupport_ < 0 )
|
---|
| 806 | convSupport_ = 3 ;
|
---|
| 807 | Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
|
---|
| 808 | convFunc.resize( convSize ) ;
|
---|
| 809 | spheroidalFunc( convFunc ) ;
|
---|
| 810 | }
|
---|
| 811 | else if ( convType_ == "GAUSS" ) {
|
---|
[2363] | 812 | // to take into account Gaussian tail
|
---|
[2356] | 813 | if ( convSupport_ < 0 )
|
---|
[2363] | 814 | convSupport_ = 12 ; // 3 * 4
|
---|
| 815 | else {
|
---|
| 816 | convSupport_ = userSupport_ * 4 ;
|
---|
| 817 | }
|
---|
[2356] | 818 | Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
|
---|
| 819 | convFunc.resize( convSize ) ;
|
---|
| 820 | gaussFunc( convFunc ) ;
|
---|
| 821 | }
|
---|
[2359] | 822 | else if ( convType_ == "PB" ) {
|
---|
| 823 | if ( convSupport_ < 0 )
|
---|
| 824 | convSupport_ = 0 ;
|
---|
[2356] | 825 | pbFunc( convFunc ) ;
|
---|
[2359] | 826 | }
|
---|
[2356] | 827 | else {
|
---|
| 828 | throw AipsError( "Unsupported convolution function" ) ;
|
---|
| 829 | }
|
---|
| 830 | }
|
---|
| 831 |
|
---|
| 832 | string STGrid::saveData( string outfile )
|
---|
| 833 | {
|
---|
[2368] | 834 | LogIO os( LogOrigin("STGrid", "saveData", WHERE) ) ;
|
---|
| 835 | double t0, t1 ;
|
---|
| 836 | t0 = mathutil::gettimeofday_sec() ;
|
---|
| 837 |
|
---|
[2360] | 838 | //Int polno = 0 ;
|
---|
[2371] | 839 | String outfile_ ;
|
---|
[2356] | 840 | if ( outfile.size() == 0 ) {
|
---|
| 841 | if ( infile_.lastchar() == '/' ) {
|
---|
| 842 | outfile_ = infile_.substr( 0, infile_.size()-1 ) ;
|
---|
| 843 | }
|
---|
| 844 | else {
|
---|
| 845 | outfile_ = infile_ ;
|
---|
| 846 | }
|
---|
| 847 | outfile_ += ".grid" ;
|
---|
| 848 | }
|
---|
| 849 | else {
|
---|
| 850 | outfile_ = outfile ;
|
---|
| 851 | }
|
---|
[2371] | 852 | Table tab ;
|
---|
| 853 | prepareTable( tab, outfile_ ) ;
|
---|
[2356] | 854 | IPosition dshape = data_.shape() ;
|
---|
[2361] | 855 | Int nrow = nx_ * ny_ * npol_ ;
|
---|
| 856 | tab.rwKeywordSet().define( "nPol", npol_ ) ;
|
---|
[2360] | 857 | tab.addRow( nrow ) ;
|
---|
[2356] | 858 | Vector<Double> cpix( 2 ) ;
|
---|
| 859 | cpix(0) = Double( nx_ - 1 ) * 0.5 ;
|
---|
| 860 | cpix(1) = Double( ny_ - 1 ) * 0.5 ;
|
---|
| 861 | Vector<Double> dir( 2 ) ;
|
---|
| 862 | ArrayColumn<Double> directionCol( tab, "DIRECTION" ) ;
|
---|
| 863 | ArrayColumn<Float> spectraCol( tab, "SPECTRA" ) ;
|
---|
| 864 | ScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
|
---|
| 865 | Int irow = 0 ;
|
---|
[2371] | 866 | Vector<Float> sp( nchan_ ) ;
|
---|
| 867 | Bool bsp, bdata ;
|
---|
| 868 | const Float *data_p = data_.getStorage( bdata ) ;
|
---|
| 869 | Float *wsp_p, *sp_p ;
|
---|
| 870 | const Float *wdata_p = data_p ;
|
---|
| 871 | long step = nx_ * ny_ * npol_ ;
|
---|
| 872 | long offset ;
|
---|
[2356] | 873 | for ( Int iy = 0 ; iy < ny_ ; iy++ ) {
|
---|
[2371] | 874 | dir(1) = center_(1) - ( cpix(1) - (Double)iy ) * celly_ ;
|
---|
[2356] | 875 | for ( Int ix = 0 ; ix < nx_ ; ix++ ) {
|
---|
[2371] | 876 | dir(0) = center_(0) - ( cpix(0) - (Double)ix ) * cellx_ ;
|
---|
[2361] | 877 | for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
|
---|
[2371] | 878 | offset = ix + iy * nx_ + ipol * nx_ * ny_ ;
|
---|
| 879 | //os << "offset = " << offset << LogIO::POST ;
|
---|
| 880 | sp_p = sp.getStorage( bsp ) ;
|
---|
| 881 | wsp_p = sp_p ;
|
---|
| 882 | wdata_p = data_p + offset ;
|
---|
| 883 | for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
|
---|
| 884 | *wsp_p = *wdata_p ;
|
---|
| 885 | wsp_p++ ;
|
---|
| 886 | wdata_p += step ;
|
---|
| 887 | }
|
---|
| 888 | sp.putStorage( sp_p, bsp ) ;
|
---|
[2356] | 889 | spectraCol.put( irow, sp ) ;
|
---|
| 890 | directionCol.put( irow, dir ) ;
|
---|
[2360] | 891 | polnoCol.put( irow, pollist_[ipol] ) ;
|
---|
[2356] | 892 | irow++ ;
|
---|
| 893 | }
|
---|
| 894 | }
|
---|
| 895 | }
|
---|
[2371] | 896 | data_.freeStorage( data_p, bdata ) ;
|
---|
[2368] | 897 |
|
---|
| 898 | t1 = mathutil::gettimeofday_sec() ;
|
---|
| 899 | os << "saveData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
|
---|
[2356] | 900 |
|
---|
| 901 | return outfile_ ;
|
---|
| 902 | }
|
---|
| 903 |
|
---|
[2371] | 904 | void STGrid::prepareTable( Table &tab, String &name )
|
---|
| 905 | {
|
---|
| 906 | Table t( infile_, Table::Old ) ;
|
---|
| 907 | t.deepCopy( name, Table::New, False, t.endianFormat(), True ) ;
|
---|
| 908 | tab = Table( name, Table::Update ) ;
|
---|
[2356] | 909 | }
|
---|
[2371] | 910 | }
|
---|