trkacxfit.cc 7.07 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 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 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 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
/*--- 
   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 
   Example command to run : 
  csh> ../AnaPaon4/Objs/trkacxfit -D 4.5 -sdec 10 -out A21_fac.txt dt_A21,trk_A21_43057,1358,1410,1278.5  dt_A21,trk_A21_43055,1512,1555,1278.5

 ----------------------------------------------------------*/

#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 

static vector<string> trksetfiles;   // datacard files defining the track/data sets 

//---   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;

    vector<AcxDataSet> v_acxd;
    vector<TrackSet> v_trk;
    for(size_t i=0; i<trksetfiles.size(); i++) {
      TrkInputDataSet tids(trksetfiles[i]);
      AcxDataSet  acxd(tids);
      TrackSet  tks(tids);
      ACxSetFitter acxfit(acxd, tks);
      string acofile = "ac_"+outfilename;
      acxfit.doACfit(acofile);  
      string acckfile = "ac_"+checkfilename;
      acxfit.saveExpectedAC(acckfile);
      string cxofile = "cx_"+outfilename;
      string cxckfile = "ac_"+checkfilename;
      acxfit.doCxfit(cxofile, cxckfile, fguseAac4Cx); 
      
      v_acxd.push_back(acxd);
      v_trk.push_back(tks);
    }
    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"
	 << "  options: -D dish_diameter -ngb -out OutBaseFileName -ckf CheckBaseFileName \n" 
	 << "           -prt PrintLevel  \n"<<endl;
    if (!fglonghelp) {
      cout << " trkacxfit -h  to get option description " << endl;
      return 2;
    }
    cout << "  -D dish_diameter : define effective dish diameter (in meter D*eff , def=4.5) \n"
      	 << "  -ngb : Use non gaussian beam profile (Bessel j1(angle))  \n"
	 << "  -prt PrintLevel: specify print level \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"<<endl;
    return 1;
  }

  outfilename = "trkfit.txt";
  checkfilename = "chktrkfit.ppf";
  prtlevel=0;
  D_dish = 4.5;
  fggaussbeam=true;

  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(); 
    }
    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;
  }

  cout << " ------------------- trkacxfit/run parameters:"<<endl; 
  cout << " Beam: D_dish (initial guess)="<<D_dish<<" GaussianBeam ? "<<(fggaussbeam?"Yes":"No")<<endl;
  cout << " OutFileName= "<<outfilename<<" CheckFileName= "<<checkfilename<<endl;
  cout << " --- TrackSetFiles: (NbFiles="<<lastargs.size()<<endl;
  trksetfiles = lastargs;
  for (size_t i=0; i<trksetfiles.size(); i++) {
    cout << "["<<i<<"] : "<< trksetfiles[i]<<endl;
  }
  cout<<"---------------------------------------------------------------------------------"<<endl;
  return 0;
}