trkacxfit.cc 13.1 KB
Newer Older
1
/*---------------------------- 
2 3 4 5 6 7 8 9
   PAON4 analysis : 
   Determining antenna pointing and dish size fitting parameters ... 
   (source track corresponds to a satellite track (or planet or sun ...)
   Does also fit of phases for the 6 cross-correlations
   Uses gacfit.h and gcxfit.h acbeam.h cxbeam.h  
   R. Ansari  , February 2019 

   Evolution of trkacfit.cc , intended to replace it 
10 11 12 13 14 15

   Assuming that the file trk_1.d contains :
@zenang 10 
@trk  dt_A21  1358,1410  1278.5  trk_A21_43057
@trk  dt_A21  1512,1555  1278.5  trk_A21_43055    

16
   Example command to run : 
17
csh> ../AnaPaon4/Objs/trkacxfit -D 4.5 -out A21_fac.txt  trk_1.d  
18

19 20 21 22 23 24 25 26 27 28

  visi-track-set datacard format 
----------------------------------
# input path (directory for files) listed here (optional)
@inpath  path 
## Zenith angle or declination shift 
# Zenith angle_in_degree (positive toward north)  initial value for fit 
@zenang ZenithAngle 
### Input data and track files - multiple @trk cards 
# filenames with .ppf extension , Tstart,Tend in minutes 
29 30
@trk AutoCrossVisiDataTableFile  Tstart,Tend  Freq  TrackFile [NOAC NOCX  NOACX]
### Optional last argument [ NOAC NOCX  NOACX ] not used currently 
31
 ----------------------------------------------------------------------------*/
32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88

#include <stdlib.h>
#include <string.h>

#include <iostream>
#include <fstream>

#include <cmath>
#include <complex>

#include "array.h"
#include "randr48.h"
#include "skymap.h"
#include "samba.h"
#include "fftwserver.h"
#include "strutilxx.h"

#include "ctimer.h"
#include "tarrinit.h"
#include "skymapinit.h"
#include "fitsioserver.h"
#include "slininterp.h"

#include "p4autils.h"

#include "acbeam.h"
#include "gacfit.h"
#include "gcxfit.h"
#include "gcxfitbaseline.h"
#include "trkfit.h"

using namespace std; 
using namespace SOPHYA; 

//---------------------------------------------------------------------------
//---  Grouping here beam & computation parameters -------------
//--- Global parameters 

// static double latitude=45.;   // latitude in degree 
// static double longitude=0.;   // longitude in degree 


static double D_dish = 4.5;       // effective dish size = effD*D 
static bool fggaussbeam=true;

// --- if true, perform Beam normalisation 
// static bool fgdobeamnormalise=false;

static string outfilename;            // output filename for auto-correlation fit parameters  
static string checkfilename;       // output filename for PAON4 and computed auto-cor using fitted parameters  

static bool fguseAac4Cx=true;       //if true, use fitted Amplitudes on autocor for initial Cross-cor fit parameter value

// static bool fgoutfits = false;   // true -> write output in fits format 

static int prtlevel=0;  // print level 

89
static string default_input_path;
90 91
static vector<string> trksetfiles;   // datacard files defining the track/data sets 

92
static bool do_cxfit = true ;       // if false, perform AutoCor fit only
93 94
// fg_cxB0: if true, fix offset B to zero (complex(0,0)) when fitting cross-correlations , false : fit B0 
static bool fg_cxB0 = true ;       // if false, fit the  offset B  when fitting cross-correlations 
95 96
static bool fg_phi0only = true;     // if false, linear varying phase with frequency, in AutoCor and Cross-Cor fits

97
// ---- fit simultane sur les 6 cross-correlations 
98
static bool do_baselinefit = false;   // if true , perform baseline fit 
99
static bool fg_fixbaseline = false;   // if true , do phases fit with fixed baselines 
100
static bool do_baselinesimplex = false;   // if true , perform baseline determination using Simplex minimisation  
101

102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121
//---   End of list of global parameters 



//-------- Declaration of functions in this file - code after the main, below 
int decode_args(int narg, char* arg[]);


//--------------------------------------------------------------
//---------  Main program -------------
//--------------------------------------------------------------
int main (int narg, char* arg[])
{
  cout << " ------ trkacxfit.cc : Fitting auto-correlations using track (satellites...) files  ------- " << endl;
  int rc = 0;
  try {
    SophyaInit();
    int rcda=decode_args(narg, arg);
    if (rcda)  return rcda;

122
    TrkFit_SetPrintLevel(prtlevel);
123
    TrkFit_FitLibInfo(); 
124 125 126
    vector<AcxDataSet> v_acxd;
    vector<TrackSet> v_trk;
    for(size_t i=0; i<trksetfiles.size(); i++) {
127 128 129
      TrkInputDataSet tids(trksetfiles[i], default_input_path);
      cout<<" --- Processing data-track set # "<<i+1<<" /max="<<trksetfiles.size()<<endl;
      cout<<tids;
130 131 132
      AcxDataSet  acxd(tids);
      TrackSet  tks(tids);
      ACxSetFitter acxfit(acxd, tks);
133 134 135
      char fnbuff[32];
      sprintf(fnbuff,"ac_%d_",(int)(i+1));
      string acofile = fnbuff+outfilename;
136
      acxfit.doACfit(acofile);  
137
      string acckfile = fnbuff+checkfilename;
138
      acxfit.saveExpectedAC(acckfile);
139 140 141 142
      if (do_cxfit)  {  //  Fit des cross-cor aussi 
	sprintf(fnbuff,"cx_%d_",(int)(i+1));
	string cxofile = fnbuff+outfilename;
	string cxckfile = fnbuff+checkfilename;
143
	acxfit.doCxfit(cxofile, fguseAac4Cx, fg_cxB0, fg_phi0only); 
144 145
	acxfit.saveExpectedCx(cxckfile);
      }
146 147 148
      v_acxd.push_back(acxd);
      v_trk.push_back(tks);
    }
149

150 151 152 153 154 155 156 157
    if (do_cxfit)  {  //  Resume resultat fit cross-cor
      cout<<" ----------  Summary phase fitting for cross-correlations ---------- " << endl;
      for(size_t i=0; i<v_acxd.size(); i++) {
	cout << "================================= TrackSet["<<i+1<<"] ================================="<<endl;
	v_acxd[i].PrintCxPhaseFitSummary(cout);
      }
      cout<<"=================================================================================="<<endl;
    }
158 159 160 161 162 163
    if (do_baselinefit || do_baselinesimplex)  {   // fit simultane des 6 cross-cor 
      string cxbofile = "cxb6_"+outfilename;
      string cxbckfile = "cxb6_"+checkfilename;
      CxBaselineFitter cxbfit(v_acxd, v_trk);
      //cxbfit.doSimplexMinimize();
      cxbfit.initFitParams();
164
      cxbfit.saveExpectedCx(cxbckfile);
165 166 167 168 169 170 171
      if (do_baselinesimplex) {
	cxbfit.doSimplexMinimize();
	cxbckfile="cxb6s_"+checkfilename;
	cxbfit.saveExpectedCx(cxbckfile);
	//      cxbfit.doCheck();
      }
      if (do_baselinefit) {
172
	cxbfit.dofit(cxbofile,fg_fixbaseline,fg_phi0only);
173 174 175
	cxbckfile = "cxb6f_"+checkfilename;
	cxbfit.saveExpectedCx(cxbckfile);
      }
176
    }
177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210
    Timer tm("trkacxfit");
  }  // End of try bloc 
  catch (PThrowable & exc) {
    cerr << " trkacxfit.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name() 
         << " - Msg= " << exc.Msg() << endl;
    rc = 99;
  }
  catch (std::exception & e) {
    cerr << " trkacxfit.cc: Catched std::exception "  
         << " - what()= " << e.what() << endl;
    rc = 98;
  }
  catch (...) {
    ;    cerr << " trkacxfit.cc: some other exception (...) was caught ! " << endl;
    rc = 97;
  }
  cout << " ------------ END OF trkacxfit.cc  (Rc= " << rc << ") ------------ " << endl;
  return rc;
}


//-----------------------------------------------------------------------------------------
//-----------------------------------------------------------------------------------------


/* --Fonction-- */
int decode_args(int narg, char* arg[])
//! Decode program arguments 
{
  bool fglonghelp=false;
  if ((narg>1)&&(strcmp(arg[1],"-h")==0))  fglonghelp=true;
  if ((narg<2)||fglonghelp) {
    cout << "  trkacxfit : fit array geometry and   \n"
	 << "  Usage: trkacxfit [-options]  track_Set_1 [track_Set_2]  [track_Set_3 ...] \n"
211
	 << "  options: -inp def_input_path -out OutFileName -ckf CheckFileName \n" 
212
	 << "           -nocxf -cxfB -phifreq -docx6f -fixb -docx6s  \n"
213
	 <<"            -D dish_diameter -ngb  -prt PrintLevel  \n"<<endl;
214 215 216 217
    if (!fglonghelp) {
      cout << " trkacxfit -h  to get option description " << endl;
      return 2;
    }
218 219 220 221 222 223
    cout << "  -inp def_input_path : default directory path for input files (track_Set_J files) \n"
      	 << "    This path can be overwritten in each track_Set_J datacard file  \n"
      	 << "  -out OutFileName : Output file (text) to save fitted parameters vectors (def trkfit.txt) \n"
      	 << "  -ckf CheckFileName : Output file (PPF) to save expected signals (def chktrkfit.ppf)  \n"
	 << "    Output file names are constructed from OutFileName and CheckFileName prepending \n" 
	 << "  ac_JJ_ or cx_JJ_  where JJ=1...n is the track_set number \n"
224
      	 << "  -nocxf : Don't perform fit on cross-cor (AutoCor only fit) - default YES -> DO Cx fits) \n"
225
      	 << "  -cxfB : Fit the offset B when fitting cross-cor (default NO -> fix B0=complex(0,0)) \n"
226 227
	 << "  -phifreq : Use linearly varying phase with frequency model - default No fixed Phi=Phi0  \n"
	 <<"             Linear phi model : Phi(freq)=Phi_0 + a_Phi*(freq-1250)/250 \n"
228 229
	 << "  -docx6f : Perform simultaneous fit over 6 cross-corr (baseline and phases determination) - default NO  \n"
	 << "  -fixb : perform the previous fit (docx6f) with fixed baselines - default NO  \n"
230
	 << "  -docx6s : try simplex minimisation over 6 cross-corr (baseline and phases determination) - default NO  \n"
231
	 << "  -ngb : Use non gaussian beam profile (Bessel j1(angle) - default Gaussian beam)  \n"
232
	 << "  -D dish_diameter : define effective dish diameter (in meter D*eff , def=4.5) \n"
233
	 << "  -prt PrintLevel: specify print level \n" 
234 235 236 237
	 << "    track_Set_J (J=1..n) are the auto/cross correlation and track data sets  \n" 
	 << "    in datacard format [ @trk (multiple) @zenang cards @inpath] \n" 
	 << "    @trk AutoCrossVisDataTableFile Tstart,Tend Freq TrackFileName \n" 
	 << "    Filenames without .ppf extension, Tstart,Tend in minutes Freq in MHz for @trk cards \n" <<endl;
238 239 240 241 242 243 244 245
    return 1;
  }

  outfilename = "trkfit.txt";
  checkfilename = "chktrkfit.ppf";
  prtlevel=0;
  D_dish = 4.5;
  fggaussbeam=true;
246 247
  do_baselinefit=false;
  fg_fixbaseline=false;
248
  do_baselinesimplex=false;
249
  fg_cxB0=true;
250 251 252 253 254 255 256 257 258 259 260

  vector<string> lastargs;
  while (narg>1) {
    string fbo = arg[1];
    if (fbo=="-D") {  // dish size   (effective dish diameter) 
      if (narg<2) { cout << "trkacxfit/decode_args missing/bad argument, -h for help " << endl;  return -1; }
      D_dish=atof(arg[2]);   arg+=2; narg-=2;   lastargs.clear(); 
    }
    else if (fbo=="-ngb") {  // Use Non gaussian beam 
      fggaussbeam=false;   arg++;  narg--;  lastargs.clear(); 
    }
261 262 263
    else if (fbo=="-nocxf") {  // DON'T FIT Cross-Correlations 
      do_cxfit=false;   arg++;  narg--;  lastargs.clear(); 
    }
264 265
    else if (fbo=="-cxfB") {  // Determine B  when fitting Cross-cor
      fg_cxB0=false;   arg++;  narg--;  lastargs.clear(); 
266
    }
267 268 269
    else if (fbo=="-phifreq") {  // Use linear model for phase(frequency)  - If not Phi(freq) = Phi0  
      fg_phi0only=false;   arg++;  narg--;  lastargs.clear(); 
    }
270
    else if (fbo=="-docx6f") {  // Perform  simultaneous fit over 6 cross-corr 
271
      do_baselinefit=true;   arg++;  narg--;  lastargs.clear(); 
272
    }
273 274 275
    else if (fbo=="-fixb") {  // Perform  the previous docx6f fit with fixed baselines (phases only)
      fg_fixbaseline=true;   arg++;  narg--;  lastargs.clear(); 
    }
276 277 278
    else if (fbo=="-docx6s") {  // Perform  simultaneous fit over 6 cross-corr 
      do_baselinesimplex=true;   arg++;  narg--;  lastargs.clear(); 
    }
279 280 281 282
    else if (fbo=="-inp") {  // output file name 
      if (narg<2) { cout << "trkacxfit/decode_args missing/bad argument, -h for help " << endl;  return -1; }
      default_input_path=arg[2];   arg+=2; narg-=2;  lastargs.clear(); 
    }
283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302
    else if (fbo=="-out") {  // output file name 
      if (narg<2) { cout << "trkacxfit/decode_args missing/bad argument, -h for help " << endl;  return -1; }
      outfilename=arg[2];   arg+=2; narg-=2;  lastargs.clear(); 
    }
    else if (fbo=="-ckf") {  // output file name for checking (data & fitted 
      if (narg<2) { cout << "trkacxfit/decode_args missing/bad argument, -h for help " << endl;  return -1; }
      checkfilename=arg[2];   arg+=2; narg-=2;  lastargs.clear(); 
    }
    else if (fbo=="-prt") {  // print level 
      if (narg<2) { cout << "trkacxfit/decode_args missing/bad argument, -h for help " << endl;  return -1; }
      prtlevel=atoi(arg[2]);   arg+=2; narg-=2;  lastargs.clear(); 
    }
    else { arg++; narg--; lastargs.push_back(fbo); }
  }

  if (lastargs.size() < 1) {
    cout << "  trkacxfit/ERROR : No input track/data specified (trkacxfit -h for usage)  \n";
    return 9;
  }

303 304
  if (!do_cxfit) do_baselinefit=do_baselinesimplex=false;
  
305 306 307
  cout << " ------------------- trkacxfit/run parameters:"<<endl; 
  cout << " Beam: D_dish (initial guess)="<<D_dish<<" GaussianBeam ? "<<(fggaussbeam?"Yes":"No")<<endl;
  cout << " OutFileName= "<<outfilename<<" CheckFileName= "<<checkfilename<<endl;
308
  cout << " Perform fits on cross-correlations ? " << (do_cxfit?"Yes":"No")
309
       << "  Phase model : Phase(freq)= " << (fg_phi0only?"Phi0":"Phi0+aPhi*(freq-1250)/250") 
310
       << "  Fit B0 ? "<<(fg_cxB0?"No (B=(0,0)":"Yes")<<endl;
311
  cout << " Perform baseline fit ? " << (do_baselinefit?"Yes":"No")
312 313
       <<(fg_fixbaseline?"  With FIXED baselines":" (Phases & baselines)")
       <<" Simplex-Minimisation ? "<<(do_baselinesimplex?"Yes":"No")<<endl;
314
  cout << " --- TrackSetFiles: (NbFiles="<<lastargs.size()<<" default directory: "<<default_input_path<<")"<<endl;
315 316 317 318 319 320 321 322
  trksetfiles = lastargs;
  for (size_t i=0; i<trksetfiles.size(); i++) {
    cout << "["<<i<<"] : "<< trksetfiles[i]<<endl;
  }
  cout<<"---------------------------------------------------------------------------------"<<endl;
  return 0;
}