Commit 85ce088e authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

Premiere version fonctionnelle du programme de fit des parametres de geometrie...

Premiere version fonctionnelle du programme de fit des parametres de geometrie PAON4 a partir des signaux Satellites (+ciel), Reza , 13/12/2018
parent 365139dc
......@@ -6,10 +6,164 @@
#define GACFIT_H_SEEN
#include "generalfit.h"
#include "simplex.h"
#include "acbeam.h"
class
static int debug_galfit=0;
//-----------------------------------------------------------------
class MyACSignal {
public:
MyACSignal(vector< vector<double> > & v_time_data_, vector< vector< vector<double> > > & vv_data_,
vector< vector< vector<double> > > & vv_err_,
vector< SLinInterp1D > & v_interp_elev_, vector< SLinInterp1D > & v_interp_cphi_,
size_t nac=0, double lambda=0.23, bool fggauss_beam=true)
: v_time_data(v_time_data_), vv_data(vv_data_), vv_err(vv_err_),
v_interp_elev(v_interp_elev_), v_interp_cphi(v_interp_cphi_),
nac_(nac), lambda_(lambda), fggauss_beam_(fggauss_beam)
{
cout << " MyACSignal() Antenna=nac+1= "<<nac_+1<<" NPar="<< 3+2*v_time_data_.size() << endl;
}
virtual int getDataSignal(size_t j, vector<double> & sig)
{
sig.resize(v_time_data[j].size());
vector<double> & vdata = vv_data[j][nac_];
for(size_t k=0; k<sig.size(); k++) sig[k]=vdata[k];
return (int)sig.size();
}
virtual TVector<double> getDataSignal(size_t j)
{
vector<double> sig;
getDataSignal(j, sig);
TVector<double> rvs(sig.size());
for(size_t k=0; k<sig.size(); k++) rvs(k)=sig[k];
return rvs;
}
virtual int getExpectedSignal(size_t j, vector<double> & sig, double D_dish, double theta_b, double phi_b, double A, double B)
{
//DBG if (debug_galfit > 0) cout<<"*DBG*A*getExpectedSignal() j="<<j<<" size="<<v_time_data[j].size()<<endl;
sig.resize(v_time_data[j].size());
vector<double> & vtime = v_time_data[j];
vector<double> & vdata = vv_data[j][nac_];
vector<double> & verr = vv_err[j][nac_];
for(size_t k=0; k<vtime.size(); k++) {
double elev=v_interp_elev[j].YInterp(vtime[k]);
double thetasrc=Angle(90.-elev,Angle::Degree).ToRadian();
double phisrc=acos(v_interp_cphi[j].YInterp(vtime[k]));
if (phisrc<0.) phisrc+=Angle::TwoPiCst();
ACBeam beam(D_dish, theta_b, phi_b, lambda_);
beam.setGaussianLobe(fggauss_beam_);
double acb=beam.Value(thetasrc, phisrc);
sig[k] = A*acb+B;
//DBG cout << " *DBG* "<<k<<" Theta,Phi-src = " << Angle(thetasrc).ToDegree() << " , " << Angle(phisrc).ToDegree()
//DBG << " elev="<<elev<<" phisrc="<<Angle(phisrc).ToDegree()<<" --> acb= "<<acb<< " sig="<<sig[k]<<endl;
}
return (int)sig.size();
}
virtual int getExpectedSignal(size_t j, vector<double> & sig, const double* parm)
{
/*
param[0] = D (Dish diameter)
parm[1] = delta_elev = theta
parm[2] = azim = phi
parm[3+2*j] = Aj (Amplitude) for track j
parm[4+2*j] = Bj (Baseline / Offset) for track j
*/
double D_dish = parm[0];
double theta_b = parm[1];
double phi_b = parm[2];
double A = parm[3+2*j];
double B = parm[4+2*j];
return getExpectedSignal(j, sig, D_dish, theta_b, phi_b, A, B);
}
virtual TVector<double> getExpectedSignal(size_t j, double D_dish, double theta_b, double phi_b, double A, double B)
{
vector<double> sig;
//DBG cout << "*DBG*getExpectedSignal() : calling getExpectedSignal(j="<<j<<" ,sig,...)"<<endl;
//DBG debug_galfit=1;
getExpectedSignal(j, sig, D_dish, theta_b, phi_b, A, B);
TVector<double> rvs(sig.size());
for(size_t k=0; k<sig.size(); k++) rvs(k)=sig[k];
return rvs;
}
virtual double getXi2(const double* parm, int & npts)
{
double xi2=0.;
npts=0;
for(size_t j=0; j<v_time_data.size(); j++) {
vector<double> sig;
npts += getExpectedSignal(j, sig, parm);
vector<double> & vdata = vv_data[j][nac_];
vector<double> & verr = vv_err[j][nac_];
for(size_t k=0; k<vdata.size(); k++) {
double xx=(sig[k]-vdata[k])/verr[k];
xi2 += (xx*xx);
}
}
return xi2;
}
vector< vector<double> > & v_time_data;
vector< vector< vector<double> > > & vv_data;
vector< vector< vector<double> > > & vv_err;
vector< SLinInterp1D > & v_interp_elev;
vector< SLinInterp1D > & v_interp_cphi;
size_t nac_;
double lambda_;
bool fggauss_beam_;
};
//-----------------------------------------------------------------
class MyACGenXi2 : public GeneralXi2, MyACSignal {
public:
MyACGenXi2(vector< vector<double> > & v_time_data_, vector< vector< vector<double> > > & vv_data_,
vector< vector< vector<double> > > & vv_err_,
vector< SLinInterp1D > & v_interp_elev_, vector< SLinInterp1D > & v_interp_cphi_,
size_t nac=0, double lambda=0.23, bool fggauss_beam=true)
: GeneralXi2(3+2*v_time_data_.size()) ,
MyACSignal(v_time_data_, vv_data_, vv_err_, v_interp_elev_, v_interp_cphi_, nac, lambda, fggauss_beam)
{
}
virtual double Value(GeneralFitData& data, double* parm, int& ndataused)
/*
param[0] = D (Dish diameter)
parm[1] = delta_elev = theta
parm[2] = azim = phi
parm[3+2*j] = Aj (Amplitude) for track j
parm[4+2*j] = Bj (Baseline / Offset) for track j
*/
{
return getXi2(parm, ndataused);
}
};
//-----------------------------------------------------------------
class MyACMinZ : public MinZFunction, MyACSignal {
public:
MyACMinZ(vector< vector<double> > & v_time_data_, vector< vector< vector<double> > > & vv_data_,
vector< vector< vector<double> > > & vv_err_,
vector< SLinInterp1D > & v_interp_elev_, vector< SLinInterp1D > & v_interp_cphi_,
size_t nac=0, double lambda=0.23, bool fggauss_beam=true)
: MinZFunction(3+2*v_time_data_.size()) ,
MyACSignal(v_time_data_, vv_data_, vv_err_, v_interp_elev_, v_interp_cphi_, nac, lambda, fggauss_beam)
{
}
virtual double Value(double const parm[])
{
int ndataused;
return getXi2(parm, ndataused);
}
};
#endif
......@@ -6,11 +6,11 @@ OBJ = ./Objs/
EXE = ./Objs/
# List of include files of this package, and .o files to handle dependencies
MYINCLISTHERE = p4autils.h visip4reader.h visp4winreader.h p4gnugain.h p4gvcor.h p4phcor.h p4freqselmgr.h acbeam.h gacfit.h
MYINCLISTHERE = p4autils.h visip4reader.h visp4winreader.h p4gnugain.h p4gvcor.h p4phcor.h p4freqselmgr.h
MYOLISTHERE = $(OBJ)/p4autils.o $(OBJ)/visip4reader.o $(OBJ)/visp4winreader.o $(OBJ)/p4gnugain.o $(OBJ)/p4gvcor.o $(OBJ)/p4phcor.o $(OBJ)/p4freqselmgr.o
# Define our target list
all : rdvisip4 visi2ntac visi2dtacx visi2tmfreq tstutls p4conv2fits msvis2dt visiavg filt_blind ckconvphases vis2ra p4vdblist rdthermrfi visamm trk_ac_fit
all : rdvisip4 visi2ntac visi2dtacx visi2tmfreq tstutls p4conv2fits msvis2dt visiavg filt_blind ckconvphases vis2ra p4vdblist rdthermrfi visamm trkacfit
clean :
rm -f $(EXE)/rdvisip4 $(EXE)/visi2ntac $(EXE)/visi2dtacx $(EXE)/visi2tmfreq $(EXE)/p4conv2fits $(EXE)/msvis2dt $(EXE)/visiavg \
......@@ -90,7 +90,7 @@ trkacfit : $(EXE)/trkacfit
$(EXE)/trkacfit : $(OBJ)/trkacfit.o $(MYOLISTHERE)
$(CXXLINK) -o $(EXE)/trkacfit $(OBJ)/trkacfit.o $(MYOLISTHERE) $(SOPHYAEXTSLBLIST)
$(OBJ)/trkacfit.o : trkacfit.cc $(MYINCLISTHERE)
$(OBJ)/trkacfit.o : trkacfit.cc $(MYINCLISTHERE) acbeam.h gacfit.h
$(CXXCOMPILE) -o $(OBJ)/trkacfit.o trkacfit.cc
......
......@@ -3,7 +3,6 @@
Determining antenna pointing and dish size through scanning parameters ...
(source track corresponds to a satellite track (or planet or sun ...)
R. Ansari , December 2018
make SGP4INC=${HERESGP4}/sgp4/libsgp4/ SGP4LIB=${HERESGP4}/libsgp4/
----------------------------------------------------------*/
#include <stdlib.h>
......@@ -29,6 +28,7 @@ make SGP4INC=${HERESGP4}/sgp4/libsgp4/ SGP4LIB=${HERESGP4}/libsgp4/
#include "slininterp.h"
#include "acbeam.h"
#include "gacfit.h"
using namespace std;
......@@ -79,70 +79,96 @@ int decode_args(int narg, char* arg[]);
//--------------------------------------------------------------
int main (int narg, char* arg[])
{
cout << " ------ trk2vis.cc : Compute visibilities from a source list ------- " << endl;
cout << " ------ trkacfit.cc : Fitting auto-correlations using track (satellites...) files ------- " << endl;
int rc = 0;
try {
SophyaInit();
size_t NTRK = 2;
string dataflnm[2] = {"dt_A21.ppf", "dt_A21.ppf"};
double tstart[2]={81500., 90750.};
double tend[2]={84500., 93250.};
string satflnm[2] = {"trk_A21_43057.ppf","trk_A21_43055.ppf"};
double zenangs[2] = {11.42,11.42};
double zenang = 11.0; // 11.42;
/*
int rcda=decode_args(narg, arg);
if (rcda) return rcda;
Timer tm("trk_ac_fit");
Timer tm("trkacfit");
*/
freq0=1278.5;
double clight = PhysQty::c().SIValue();
lambda0 = clight/(freq0*1.e6);
vector< vector<double> > v_time_data(2);
vector< vector< vector<double> > > vv_data(2);
vector< vector<double> > v_time_data(NTRK);
vector< vector< vector<double> > > vv_data(NTRK);
vector< vector< vector<double> > > vv_err(NTRK);
vector< vector<double> > v_min_data(NTRK);
vector< vector<double> > v_max_data(NTRK);
int tot_npoints = 0; // total number of points for fit
{
string dataflnm = "dt_A21.ppf";
cout << "1- Extracting data from data file DataTable: " << dataflnm << endl;
double tstart[2]={81500., 90750.};
double tend[2]={84500., 93250.};
DataTable dt_data;
PInPersist pin(dataflnm);
pin >> dt_data;
dt_data.SetShowMinMaxFlag(true);
size_t ktime = dt_data.IndexNom("timesec");
vector<double> vtm;
dt_data.GetColumn(ktime, vtm);
const char * acname[4]={"V11","V22","V33","V44"};
vector< vector<double> > v_vac(4);
for(size_t ii=0; ii<4; ii++) { // 4 auto-correlations
size_t kac = dt_data.IndexNom(acname[ii]);
dt_data.GetColumn(kac, v_vac[ii]);
vector<double> vtmp1, vtmp2;
vv_data[0].push_back(vtmp1);
vv_data[1].push_back(vtmp1);
}
for(size_t j=0; j<2; j++) {
vector< vector<double> > & rvac = vv_data[j];
cout << "1."<<j+1<<" : For time interval (Sat"<<j+1<<") "<<tstart[j]<<" < t < "<<tend[j]<<endl;
for(size_t jj=0; jj<vtm.size(); jj++) {
if ((vtm[jj]<tstart[j])||(vtm[jj]>tend[j])) continue;
v_time_data[j].push_back(vtm[jj]);
for(size_t j=0; j<NTRK; j++) {
cout << "1."<<j+1<<" Extracting data from data file DataTable: " << dataflnm[j]<<endl
<< " ... For time interval (Trk"<<j+1<<") "<<tstart[j]<<" < t < "<<tend[j]<<endl;
DataTable dt_data;
PInPersist pin(dataflnm[j]);
pin >> dt_data;
dt_data.SetShowMinMaxFlag(true);
size_t ktime = dt_data.IndexNom("timesec");
vector<double> vtm;
dt_data.GetColumn(ktime, vtm);
vector< vector<double> > v_vac(4);
for(size_t ii=0; ii<4; ii++) { // 4 auto-correlations
size_t kac = dt_data.IndexNom(acname[ii]);
dt_data.GetColumn(kac, v_vac[ii]);
vector<double> vtmp, vetmp;
vv_data[j].push_back(vtmp);
vv_err[j].push_back(vetmp);
v_min_data[j].push_back(9.e19);
v_max_data[j].push_back(-9.e19);
}
vector< vector<double> > & v_data = vv_data[j];
vector< vector<double> > & v_err = vv_err[j];
for(size_t k=0; k<vtm.size(); k++) {
if ((vtm[k]<tstart[j])||(vtm[k]>tend[j])) continue;
v_time_data[j].push_back(vtm[k]);
for(size_t ii=0; ii<4; ii++) {
vector<double> & vac = v_vac[ii];
rvac[ii].push_back(vac[jj]);
v_data[ii].push_back(vac[k]);
v_err[ii].push_back(0.1*sqrt(fabs(vac[k]))); // calcul d'erreur, a affiner
if (vac[k]<v_min_data[j][ii]) v_min_data[j][ii]=vac[k];
if (vac[k]>v_max_data[j][ii]) v_max_data[j][ii]=vac[k];
}
}
cout << " ... Done for timeInterval " << j+1 << " data size="<<v_time_data[j].size()<<endl;
tot_npoints += v_time_data[j].size(); // total number of points for fit
cout << " ... Done for " << j+1 << " data size="<<v_time_data[j].size()<<endl;
cout << " Data-AutoCor Min[A1...A4]=";
for(size_t ii=0; ii<4; ii++) cout<<setw(10)<<v_min_data[j][ii]<< " , "; cout << endl;
cout << " Data-AutoCor Max[A1...A4]=";
for(size_t ii=0; ii<4; ii++) cout<<setw(10)<<v_max_data[j][ii]<< " , "; cout << endl;
}
}
string satflnm[2] = {"trk_A21_43057.ppf","trk_A21_43055.ppf"};
vector< vector<double> > v_time_sat(2);
vector< vector<double> > v_sat_elev(2);
vector< vector<double> > v_sat_azim(2);
vector< vector<double> > v_time_sat(NTRK);
vector< vector<double> > v_sat_elev(NTRK);
vector< vector<double> > v_sat_azim(NTRK);
vector< SLinInterp1D > v_interp_elev(NTRK);
vector< SLinInterp1D > v_interp_cphi(NTRK);
vector< SLinInterp1D > v_interp_elev(2);
vector< SLinInterp1D > v_interp_azim(2);
vector<double> v_Ddish(4);
vector<double> v_thetaant(4), v_phiant(4);
vector< vector<double> > v_A(4), v_B(4);
for(size_t ii=0; ii<4; ii++) {
v_A[ii].resize(NTRK); v_B[ii].resize(NTRK);
}
cout << "2- Extracting data from two satellite track DataTables: " << satflnm[0] << " , " << satflnm[1] << endl;
for(size_t j=0; j<2; j++) {
for(size_t j=0; j<NTRK; j++) {
cout << "2."<<j+1<<" Extracting data from satellite track DataTables: " << satflnm[j] << endl;
DataTable dt_sat;
PInPersist pin(satflnm[j]);
pin >> dt_sat;
......@@ -155,11 +181,119 @@ int main (int narg, char* arg[])
dt_sat.GetColumn(kazim, v_sat_azim[j]);
cout << "2."<<j+1<<" Done for satellite from file "<<satflnm[j]<<" -> size()="<<v_time_sat[j].size()<<endl;
v_interp_elev[j].DefinePoints(v_time_sat[j], v_sat_elev[j]);
v_interp_azim[j].DefinePoints(v_time_sat[j], v_sat_azim[j]);
vector<double> cazim(v_sat_azim[j].size());
for(size_t k=0; k<v_sat_azim[j].size(); k++) {
double phisrcdeg=90.-v_sat_azim[j][k];
if (phisrcdeg<0.) phisrcdeg+=360.;
double phisrc=Angle(phisrcdeg,Angle::Degree).ToRadian();
cazim[k]=cos(phisrc);
}
// v_interp_azim[j].DefinePoints(v_time_sat[j], v_sat_azim[j]);
v_interp_cphi[j].DefinePoints(v_time_sat[j], cazim);
cout<<" DONE Creation SLinInterp1D for elevation / azimuth ..."<<endl;
cout << v_interp_elev[j];
cout << v_interp_azim[j];
cout << v_interp_cphi[j];
}
for(size_t ii=0; ii<4; ii++) {
cout << "3."<<ii+1<<" Creating General Fit for AutoCor Antenna= " << ii+1 << endl;
GeneralFitData gdata(1, tot_npoints);
for(size_t j=0; j<NTRK; j++) {
vector< vector<double> > & v_data = vv_data[j];
vector< vector<double> > & v_err = vv_err[j];
for(size_t k=0; k<v_time_data[j].size(); k++) {
gdata.AddData1(v_time_data[j][k],v_data[ii][k],v_err[ii][k]); // Fill x, y and error on y
}
}
MyACGenXi2 gxi2( v_time_data, vv_data, vv_err, v_interp_elev, v_interp_cphi, ii, lambda0, true);
GeneralFit mFit(&gxi2);
mFit.SetData(&gdata); // connect data to the fitter , here the data is unused - gxi2 includes its data
mFit.SetMaxStep(1000);
// SetParam(int n,double value, double step,double min=1., double max=-1.);
mFit.SetParam(0,"D_dish",D_dish,0.1,D_dish*0.8,D_dish*1.2);
// mFit.SetFix(0, D_dish);
double thetaAntenne=0., phiAntenne=0.;
if (fabs(zenang)>1.e-6) {
if (zenang<0) {
thetaAntenne=Angle(-zenang,Angle::Degree).ToRadian();
phiAntenne=Angle(270.,Angle::Degree).ToRadian();
}
else {
thetaAntenne=Angle(zenang,Angle::Degree).ToRadian();
phiAntenne=Angle(90.,Angle::Degree).ToRadian();
}
}
mFit.SetParam(1,"ThetaAntenne",thetaAntenne,M_PI/1440,0.,thetaAntenne+M_PI/36.);
mFit.SetParam(2,"PhiAntenne",phiAntenne,M_PI/180.,0.,2.*M_PI);
// mFit.SetFix(1, thetaAntenne);
// mFit.SetFix(2, phiAntenne);
for(size_t j=0; j<NTRK; j++) {
char pname[32];
sprintf(pname,"A%d",(int)(j+1));
double A = v_max_data[j][ii];
mFit.SetParam(2*j+3,pname,A,A/10.,A/20,A*5);
sprintf(pname,"B%d",(int)(j+1));
double B = v_min_data[j][ii];
mFit.SetParam(2*j+4,pname,B,B/10.,B/20,B*5);
mFit.SetFix(2*j+4, B);
}
cout << "4."<<ii+1<<" Performing the fit for AutoCor Antenna= " << ii+1 << endl;
int rcfit = mFit.Fit();
if (prtlevel>1) mFit.PrintFit();
if(rcfit>0) {
cout<< "-------------------------- Result for Antenna No " << ii+1 << endl;
cout<< "---Fit resultReduce_Chisquare = " << mFit.GetChi2Red()<< " nstep="<<mFit.GetNStep() << " rc="<<rcfit<<endl;
double Dfit=mFit.GetParm(0); double err_Dfit=mFit.GetParmErr(0);
cout <<setw(16)<<"DishDiameter= "<<setw(12)<<Dfit<<" +/- "<<setw(12)<<err_Dfit<<" m."<<endl;
v_Ddish[ii]=Dfit;
double thetaant=mFit.GetParm(1); double err_thetaant=mFit.GetParmErr(1);
v_thetaant[ii]=thetaant;
cout <<setw(16)<<"ThetaAntenne= "<<setw(12)<<Angle(thetaant).ToDegree()<< " +/- "
<<setw(12)<<Angle(err_thetaant).ToDegree()<<" deg"<<endl;
double phiant=mFit.GetParm(2); double err_phiant=mFit.GetParmErr(2);
v_phiant[ii]=phiant;
cout <<setw(16)<<"PhiAntenne= "<<setw(12)<<Angle(phiant).ToDegree()<< " +/- "
<<setw(12)<<Angle(err_phiant).ToDegree()<<" deg"<<endl;
for(size_t j=0; j<NTRK; j++) {
double A=mFit.GetParm(3+2*j); double err_A=mFit.GetParmErr(3+2*j);
double B=mFit.GetParm(4+2*j); double err_B=mFit.GetParmErr(4+2*j);
cout << " Trk/Sat["<<j<<"] -> A= "<<A<<" +/- "<<err_A<<" B= "<<B<<" +/- "<<err_B<<endl;
v_A[ii][j]=A; v_B[ii][j]=B;
}
}
else {
cout << "--- Fit_Error, rc = " << rcfit << " nstep="<<mFit.GetNStep()<<endl;
if (prtlevel>0) mFit.PrintFitErr(rcfit);
}
}
POutPersist pos("trkac.ppf");
for(size_t ii=0; ii<4; ii++) {
cout << " 5. Computing DataSignal & Expected Signal for fitted params and dish "<<ii+1<<endl;
MyACSignal macs(v_time_data, vv_data, vv_err, v_interp_elev, v_interp_cphi, ii, lambda0, true);
double Ddishfit=v_Ddish[ii];
double thetafit=v_thetaant[ii];
double phifit=v_phiant[ii];
char oname[32];
for(size_t j=0; j<NTRK; j++) {
double A = v_A[ii][j];
double B = v_B[ii][j];
Vector signal = macs.getDataSignal(j);
sprintf(oname,"data_%d_%d",(int)ii+1,(int)j+1);
pos << PPFNameTag(oname)<<signal;
Vector expsignal = macs.getExpectedSignal(j, Ddishfit, thetafit, phifit, A, B);
sprintf(oname,"simu_%d_%d",(int)ii+1,(int)j+1);
pos << PPFNameTag(oname)<<expsignal;
}
}
/*
cout << " ------- Fitting AutoCorrelations (Dgeom="<<Dgeom<<" Eff="<<effD<<" ->D_dish=D*eff="<<D_dish
cout << " for lambda = "<<lambda0 << "m. Freq="<<freq0<<" MHz"<<endl ;
......@@ -286,20 +420,20 @@ int main (int narg, char* arg[])
*/
} // End of try bloc
catch (PThrowable & exc) {
cerr << " trk_ac_fit.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name()
cerr << " trkacfit.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name()
<< " - Msg= " << exc.Msg() << endl;
rc = 99;
}
catch (std::exception & e) {
cerr << " trk_ac_fit.cc: Catched std::exception "
cerr << " trkacfit.cc: Catched std::exception "
<< " - what()= " << e.what() << endl;
rc = 98;
}
catch (...) {
; cerr << " trk_ac_fit.cc: some other exception (...) was caught ! " << endl;
; cerr << " trkacfit.cc: some other exception (...) was caught ! " << endl;
rc = 97;
}
cout << " ------------ END OF trk_ac_fit.cc (Rc= " << rc << ") ------------ " << endl;
cout << " ------------ END OF trkacfit.cc (Rc= " << rc << ") ------------ " << endl;
return rc;
}
......@@ -332,8 +466,8 @@ int decode_args(int narg, char* arg[])
//! Decode program arguments
{
if ((narg<2)||((narg>1)&&(strcmp(arg[1],"-h")==0))) {
cout << " trk_ac_fit : computes visibility arrays from a list of sources \n"
<< " Usage: trk_ac_fit [-options] [dish1_pos dish2_pos dish3_pos ...] \n"
cout << " trkacfit : computes visibility arrays from a list of sources \n"
<< " Usage: trkacfit [-options] [dish1_pos dish2_pos dish3_pos ...] \n"
<< " options: [-paon4] [-lat latitude_degree] [-D dish_diameter] \n"
<< " [-sdec angle_degree] [-freq f0_MHz] [-lam lambda0_meter] \n"
<< " [-in InputTrackPPFFile] [-out OutVisi_FileName] [-prt PrintLevel] \n"<<endl;
......@@ -368,35 +502,35 @@ int decode_args(int narg, char* arg[])
setup_paon4(); arg++; narg--; lastargs.clear();
}
else if (fbo=="-lat") { // latitude
if (narg<2) { cout << "trk_ac_fit/decode_args missing/bad argument, -h for help " << endl; return -1; }
if (narg<2) { cout << "trkacfit/decode_args missing/bad argument, -h for help " << endl; return -1; }
latitude=atof(arg[2]); arg+=2; narg-=2; lastargs.clear();
}
else if (fbo=="-D") { // dish size
if (narg<2) { cout << "trk_ac_fit/decode_args missing/bad argument, -h for help " << endl; return -1; }
if (narg<2) { cout << "trkacfit/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=="-freq") { // frequency range (min,max,step) in MHz
if (narg<2) { cout << "trk_ac_fit/decode_args missing/bad argument, -h for help " << endl; return -1; }
if (narg<2) { cout << "trkacfit/decode_args missing/bad argument, -h for help " << endl; return -1; }
freq0=atof(arg[2]); lambda0 = clight/(freq0*1.e6); arg+=2; narg-=2; lastargs.clear();
}
else if (fbo=="-lam") { // frequency range (min,max,step) in MHz
if (narg<2) { cout << "trk_ac_fit/decode_args missing/bad argument, -h for help " << endl; return -1; }
if (narg<2) { cout << "trkacfit/decode_args missing/bad argument, -h for help " << endl; return -1; }
lambda0=atof(arg[2]); freq0 = clight/(lambda0)/1.e6; arg+=2; narg-=2; lastargs.clear();
}
else if (fbo=="-sdec") { // pointing declination
if (narg<2) { cout << "trk_ac_fit/decode_args missing/bad argument, -h for help " << endl; return -1; }
if (narg<2) { cout << "trkacfit/decode_args missing/bad argument, -h for help " << endl; return -1; }
shift_dec=atof(arg[2]); arg+=2; narg-=2; lastargs.clear();
}
else if (fbo=="-in") { // input track DataTable file name
if (narg<2) { cout << "trk_ac_fit/decode_args missing/bad argument, -h for help " << endl; return -1; }
if (narg<2) { cout << "trkacfit/decode_args missing/bad argument, -h for help " << endl; return -1; }
infilename=arg[2]; arg+=2; narg-=2; lastargs.clear();
}
else if (fbo=="-out") { // output file name
if (narg<2) { cout << "trk_ac_fit/decode_args missing/bad argument, -h for help " << endl; return -1; }
if (narg<2) { cout << "trkacfit/decode_args missing/bad argument, -h for help " << endl; return -1; }
outfilename=arg[2]; arg+=2; narg-=2; lastargs.clear();
}
else if (fbo=="-prt") { // print level
if (narg<2) { cout << "trk_ac_fit/decode_args missing/bad argument, -h for help " << endl; return -1; }
if (narg<2) { cout << "trkacfit/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); }
......@@ -411,11 +545,11 @@ int decode_args(int narg, char* arg[])
}
}
if (dish_pos.size() < 2) {
cout << "trk_ac_fit/decode_args()/ERROR specify at least two dish positions ... " << endl;
cout << "trkacfit/decode_args()/ERROR specify at least two dish positions ... " << endl;
return 2;
}
cout << " ---- trk_ac_fit/run parameters:"<<endl;
cout << " ---- trkacfit/run parameters:"<<endl;
cout << " Latitude="<<latitude<<" deg D_dish="<<D_dish<<" m. Frequency: "<<freq0<<" MHz "<<" -> Lambda "<<lambda0<<" m."<<endl;
cout << " Computing baselines from "<< dish_pos.size() << " dish positions"<<endl;
for (size_t i=0; i<dish_pos.size(); i++) {
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment