Source
xxxxxxxxxx
7232
7232
}
7233
7233
else
7234
7234
logSink() << " INSUFFICIENT DATA ";
7235
7235
7236
7236
logSink() << LogIO::POST;
7237
7237
*/
7238
7238
} // ispw
7239
7239
}
7240
7240
} // iTran
7241
7241
// max 3 coefficients
7242
-
//Matrix<Double> spidx(nTran,3,0.0);
7243
-
//Matrix<Double> spidxerr(nTran,3,0.0);
7244
-
Matrix<Double> spidx(nFld,3,0.0);
7245
-
Matrix<Double> spidxerr(nFld,3,0.0);
7242
+
//Matrix<Double> spidx(nFld,3,0.0);
7243
+
//Matrix<Double> spidxerr(nFld,3,0.0);
7244
+
Matrix<Double> spidx(nFld,fitorder+1,0.0);
7245
+
Matrix<Double> spidxerr(nFld,fitorder+1,0.0);
7246
7246
Matrix<Double> covar;
7247
7247
Vector<Double> fitFluxD(nFld,0.0);
7248
7248
Vector<Double> fitFluxDErr(nFld,0.0);
7249
7249
Vector<Double> fitRefFreq(nFld,0.0);
7250
7250
//Vector<Double> refFreq(nFld,0.0);
7251
7251
7252
7252
for (Int iTran=0; iTran<nTran; iTran++) {
7253
7253
uInt tranidx=tranField(iTran);
7254
7254
Int nValidFlux=ntrue(scaleOK.column(tranidx));
7255
7255
7267
7267
// shift mean(log(freq)) later
7268
7268
// Reference frequency is first in the list
7269
7269
//refFreq(tranidx)=freqs[0];
7270
7270
//freqs/=refFreq(tranidx);
7271
7271
//
7272
7272
// calculate spectral index
7273
7273
// fit the per-spw fluxes to get spectral index
7274
7274
LinearFit<Double> fitter;
7275
7275
uInt myfitorder;
7276
7276
if (nValidFlux > 2) {
7277
-
if (fitorder > 2) {
7278
-
logSink() << LogIO::WARN << "Currently only support fitorder < 3, using fitorder=2 instead"
7279
-
<< LogIO::POST;
7280
-
myfitorder = 2;
7277
+
//if (fitorder > 2) {
7278
+
// logSink() << LogIO::WARN << "Currently only support fitorder < 3, using fitorder=2 instead"
7279
+
// << LogIO::POST;
7280
+
// myfitorder = 2;
7281
+
//}
7282
+
//else {
7283
+
// if (fitorder < 0) {
7284
+
// logSink() << LogIO::WARN
7285
+
// << "fitorder=" << fitorder
7286
+
// << " not supported. Using fitorder=1"
7287
+
// << LogIO::POST;
7288
+
// myfitorder = 1;
7289
+
// }
7290
+
// else {
7291
+
// myfitorder = (uInt)fitorder;
7292
+
// }
7293
+
7294
+
if (fitorder < 0) {
7295
+
logSink() << LogIO::WARN
7296
+
<< "fitorder=" << fitorder
7297
+
<< " not supported. Using fitorder=1"
7298
+
<< LogIO::POST;
7299
+
myfitorder = 1;
7281
7300
}
7282
7301
else {
7283
-
if (fitorder < 0) {
7284
-
logSink() << LogIO::WARN
7285
-
<< "fitorder=" << fitorder
7286
-
<< " not supported. Using fitorder=1"
7287
-
<< LogIO::POST;
7288
-
myfitorder = 1;
7289
-
}
7290
-
else {
7291
-
myfitorder = (uInt)fitorder;
7292
-
}
7302
+
if (fitorder < nValidFlux) {
7303
+
myfitorder = (uInt)fitorder;
7304
+
}
7305
+
else {
7306
+
myfitorder = (uInt)(nValidFlux-1);
7307
+
logSink() << LogIO::WARN
7308
+
<< "Not enough number of valid flux density data for the request fitorder:"<<fitorder
7309
+
<<". Using a lower fitorder="<<myfitorder<<LogIO::POST;
7310
+
}
7293
7311
}
7294
7312
}
7295
7313
else {
7296
7314
myfitorder = 1;
7297
7315
}
7298
7316
// set fitting for spectral index, alpha and beta(curvature)
7299
7317
// with S = S_o*f/f0^(alpha+beta*log(f/fo))
7300
7318
// fitting y = c + alpha*x + beta*x^2
7301
7319
// log(S/S0)=alpha*log(f/f0) + beta*log(f/f0)**2
7302
7320
Polynomial< AutoDiff<Double> > bp(myfitorder);
7325
7343
oFitMsg += fldNames(tranidx);
7326
7344
oFitMsg += " with fitorder="+String::toString<Int>(myfitorder)+": ";
7327
7345
// oFitMsg += "Flux density = "+String::toString<Double>(pow(10.0,(soln(0))));
7328
7346
oFitMsg += "Flux density = "+String::toString<Double>(fitFluxD(tranidx));
7329
7347
// Double ferr=(errs(0)>0.0 ? exp10(errs(0)) : 0.0);
7330
7348
// Double ferr=(errs(0)>0.0 ? pow(10.0,(errs(0))) : 0.0);
7331
7349
oFitMsg += " +/- "+String::toString<Double>(fitFluxDErr(tranidx));
7332
7350
// oFitMsg += " (freq="+String::toString<Double>(refFreq(tranidx)/1.0e9)+" GHz)";
7333
7351
// oFitMsg += " (freq="+String::toString<Double>(pow(10.0,meanLogFreq)/1.0e9)+" GHz)";
7334
7352
oFitMsg += " (freq="+String::toString<Double>(fitRefFreq(tranidx)/1.0e9)+" GHz)";
7353
+
//oFitMsg += " soln.nelements="+String::toString<Int>(soln.nelements());
7354
+
oFitMsg += " spidx:";
7335
7355
for (uInt j=1; j<soln.nelements();j++) {
7336
-
if (j==1) {
7337
-
oFitMsg += " spidx="+String::toString<Double>(soln(1));
7338
-
if (nValidFlux>2)
7339
-
oFitMsg += " +/- "+String::toString<Double>(errs(1));
7340
-
else
7341
-
oFitMsg += " (degenerate)";
7356
+
String coefname=" a_"+String::toString<Int>(j);
7357
+
if (j==1) coefname += " (spectral index) ";
7358
+
oFitMsg += coefname+"="+String::toString<Double>(soln(j));
7359
+
if (nValidFlux > (Int)(j+1)) {
7360
+
oFitMsg += " +/- "+String::toString<Double>(errs(j));
7342
7361
}
7343
-
if (j==2) {
7344
-
oFitMsg += " curv="+String::toString<Double>(soln(2));
7345
-
if (nValidFlux>3)
7346
-
oFitMsg += " +/- "+String::toString<Double>(errs(2));
7347
-
else
7362
+
else {
7348
7363
oFitMsg += " (degenerate)";
7349
7364
}
7350
7365
}
7351
7366
if ( oListFile != "" ) {
7352
7367
ofstream oStream;
7353
7368
oStream.open( oListFile.chars(), ios::out|ios::app );
7354
7369
oStream << "#" << oFitMsg << endl << flush;
7355
7370
oStream.close();
7356
7371
}
7357
7372
logSink() << oFitMsg << LogIO::POST;
7493
7508
for (Int iFld=0; iFld<nFld; iFld++) {
7494
7509
if (MGOK[iFld]!=NULL) {
7495
7510
delete MGOK[iFld];
7496
7511
delete MG[iFld];
7497
7512
delete MG2[iFld];
7498
7513
delete MGWT[iFld];
7499
7514
delete MGVAR[iFld];
7500
7515
delete MGN[iFld];
7501
7516
delete MGNALL[iFld];
7502
7517
}
7503
-
}
7518
+
}
7504
7519
7505
7520
// cout << "done." << endl;
7506
7521
7507
7522
// Avoid this since problem with <msname>/SELECTED_TABLE (06Feb02 gmoellen)
7508
7523
/*
7509
7524
{
7510
7525
7511
7526
MeasurementSet ms(vs_->msName().before("/SELECTED"));
7512
7527
Table historytab = Table(ms.historyTableName(),
7513
7528
TableLock(TableLock::UserNoReadLocking),