Commit ca42d02f authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

Merge par Reza suite ajout programme fit autocor ds Satellites, Reza

Merge branch 'master' of gitlab.in2p3.fr:baoradio/AnaPAON4
parents cef458f8 4dce16a0
......@@ -12,7 +12,7 @@ SGP4LIB = ${EXTLIBDIR}/Satellites/sgp4/libsgp4/
SGP4CCFLG =
# Define our target list
all : predictsatsgp4 trk2dt
all : predictsatsgp4 trk2dt trk_ac_fit
clean :
rm -f $(EXE)/predictsatsgp4
......@@ -25,6 +25,9 @@ $(OBJ)/predictsatsgp4.o : predictsatsgp4.cc $(SGP4INC)/Tle.h
$(OBJ)/trk2dt.o : trk2dt.cc
$(CXXCOMPILE) -o $(OBJ)/trk2dt.o trk2dt.cc
$(OBJ)/trk_ac_fit.o : trk_ac_fit.cc
$(CXXCOMPILE) -o $(OBJ)/trk_ac_fit.o trk_ac_fit.cc
###### Compilation et link des executables
predictsatsgp4 : $(EXE)/predictsatsgp4
......@@ -34,9 +37,17 @@ $(EXE)/predictsatsgp4 : $(OBJ)/predictsatsgp4.o
$(CXXLINK) -o $(EXE)/predictsatsgp4 $(OBJ)/predictsatsgp4.o -L$(SGP4LIB) -lsgp4 -lstdc++ -lc -lm
###### Compilation et link des executables
trk2dt : $(EXE)/trk2dt
echo 'trk2dt --- made'
$(EXE)/trk2dt : $(OBJ)/trk2dt.o
$(CXXLINK) -o $(EXE)/trk2dt $(OBJ)/trk2dt.o $(SOPHYAEXTSLBLIST)
###### Compilation et link des executables
trk_ac_fit : $(EXE)/trk_ac_fit
echo 'trk_ac_fit --- made'
$(EXE)/trk_ac_fit : $(OBJ)/trk_ac_fit.o
$(CXXLINK) -o $(EXE)/trk_ac_fit $(OBJ)/trk_ac_fit.o $(SOPHYAEXTSLBLIST)
/* PAON4 analysis software
Class representing an autocorrelation beam (adapted from JSkyMap code)
R. Ansari, December 2018 */
#ifndef GACFIT_H_SEEN
#define GACFIT_H_SEEN
#include "generalfit.h"
#include "acbeam.h"
class
#endif
......@@ -277,6 +277,9 @@ int main(int narg, char *arg[])
cout <<" velocity: |V|="<<velsat<<" km/s, angular |Vangle|="<<vangmod<<", Vaz="<<vangaz<<", Valt="<<vangalt<<" deg/min"<<endl;
cout <<" closest-los: sep="<<sepmin<<" "<<datemin<<" "<<topomin<<" "<<vsatfilename[isat]<<endl;
if(sepmin>sep) cout<<"WARNING: sep="<<sep<<"< sepmin="<<sepmin<<endl;
if (sepmin<2.) cout << "GOOD SAT " <<satname << " " << vsatfilename[isat]<<" " << " closest-los: sep="<<sepmin<<" "<<datemin<<" "<<topomin << endl;
//
} catch (SatelliteException& e) {
if(lp>0) cerr<<"Sat="<<satname<<" skipped! "<<e.what()<<" "<<vsatfilename[isat]<<endl;
......
......@@ -88,7 +88,7 @@ int main(int narg, const char* arg[])
int nok = 0;
DataTableRow row = dt.EmptyRow();
ifstream is(infile);
ifstream is(infile.c_str());
string sline;
while (!is.eof()) {
sline = "";
......
/*---
PAON4 analysis :
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>
#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 "acbeam.h"
using namespace std;
using namespace SOPHYA;
//---------------------------------------------------------------------------
//--- Grouping here beam & computation parameters -------------
//--- Global parameters
static int mHeal=512;
static int Lmax=1440;
// static int mHeal=128;
// static int Lmax=443;
static double latitude=45.; // latitude in degree
static double longitude=0.; // longitude in degree
static double Dgeom = 5.; // geometric Dish size in meters
static double effD = 0.9; // dish efficiency , filling factor
static double D_dish = 4.5; // effective dish size = effD*D
static bool fggaussbeam=true;
static vector< Vector3d > dish_pos; // dish positions
static vector< Vector3d > baselines; // baselines
static vector < pair<int, int> > baselines_pair; // pair of antenna making up the baseline
static double shift_dec=0.; // shift_dec in degree (shift in pointing angle in the NS plane, positive toward noth)
// --- if true, perform Beam normalisation in the same way as in configs.cc in CreateBeamLMList() function
static bool fgdobeamnormalise=true;
static double lambda0=0.21; // wavelength, in meters
static double freq0=1420.4; // Corresponding frequency, in MHz
static string infilename; // input track PPF file name
static string outfilename; // output visibility DataTable filename
static bool fgoutfits = false; // true -> write output in fits format
static int prtlevel=0; // print level
//--- End of list of global parameters
//-------- Declaration of utility functions - code after the main, below
void setup_paon4();
int decode_args(int narg, char* arg[]);
//--------------------------------------------------------------
//--------- Main program -------------
//--------------------------------------------------------------
int main (int narg, char* arg[])
{
cout << " ------ trk2vis.cc : Compute visibilities from a source list ------- " << endl;
int rc = 0;
try {
SophyaInit();
/*
int rcda=decode_args(narg, arg);
if (rcda) return rcda;
Timer tm("trk_ac_fit");
*/
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);
{
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 ii=0; ii<4; ii++) {
vector<double> & vac = v_vac[ii];
rvac[ii].push_back(vac[jj]);
}
}
cout << " ... Done for timeInterval " << j+1 << " data size="<<v_time_data[j].size()<<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< SLinInterp1D > v_interp_elev(2);
vector< SLinInterp1D > v_interp_azim(2);
cout << "2- Extracting data from two satellite track DataTables: " << satflnm[0] << " , " << satflnm[1] << endl;
for(size_t j=0; j<2; j++) {
DataTable dt_sat;
PInPersist pin(satflnm[j]);
pin >> dt_sat;
dt_sat.SetShowMinMaxFlag(true);
size_t ktime = dt_sat.IndexNom("timesec");
dt_sat.GetColumn(ktime, v_time_sat[j]);
size_t kelev = dt_sat.IndexNom("elevation");
dt_sat.GetColumn(kelev, v_sat_elev[j]);
size_t kazim = dt_sat.IndexNom("azimuth");
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]);
cout<<" DONE Creation SLinInterp1D for elevation / azimuth ..."<<endl;
cout << v_interp_elev[j];
cout << v_interp_azim[j];
}
/*
cout << " ------- Fitting AutoCorrelations (Dgeom="<<Dgeom<<" Eff="<<effD<<" ->D_dish=D*eff="<<D_dish
cout << " for lambda = "<<lambda0 << "m. Freq="<<freq0<<" MHz"<<endl ;
if (fggaussbeam) cout << " Using Gaussian beam ..."<<endl;
if (fgrectlobe) cout << " Using rectangular dish Lx="<<Lx<<" Ly="<<Ly<< endl;
DataTable dtk;
cout << "1- Reading input track file "<<infilename<<endl;
{
PInPersist pis(infilename);
pis >> dtk;
}
dtk.SetShowMinMaxFlag(true);
cout << dtk ;
vector< TimeStamp > vdate; dtk.GetColumn(dtk.IndexNom("datetime"), vdate);
vector< r_8 > vtms; dtk.GetColumn(dtk.IndexNom("timesec"), vtms);
vector< r_8 > vazim; dtk.GetColumn(dtk.IndexNom("azimuth"), vazim);
vector< r_8 > velev; dtk.GetColumn(dtk.IndexNom("elevation"), velev);
cout << "2- Computing beam renormalisation factor ..."<<endl;
// Vector3d DeltaR0=coortr.Transform(0.,0.,0.); // the computation is independent of the baseline and the pointing
//---- Systeme de coordonnees local Ox vers l'Est, Oy vers le Nord , Oz vertical, vers le zenith
// Angles dans le repere local , Elevation (par rapport au plan xOy) Azimuth (par rapport l'axe S-N) ,
// dans le sens croissant N->E->S-W
// Pointage de l'antenne dans le plan N-S (en declinaison)
// Si shift_dec>0 , azimuth=0 -> phi_antenne=90 , theta=shift_dec
// shift_dec<0 , azimuth=180 -> phi_antenne=270 (ou -90) , theta=-shift_dec
// Position angulaire des sources :
// - angle theta (avec Oz) = 90.-Elevation-source (degres)
// - angle phi (avec Ox) = 90.-Azimuth (degres)
double thetaAntenne = 0.;
double phiAntenne = 0.;
if (fabs(shift_dec)>1.e-6) { // non zero pointing shift
if (shift_dec<0) {
thetaAntenne=Angle(-shift_dec,Angle::Degree).ToRadian();
phiAntenne=Angle(270.,Angle::Degree).ToRadian();
}
else {
thetaAntenne=Angle(shift_dec,Angle::Degree).ToRadian();
phiAntenne=Angle(90.,Angle::Degree).ToRadian();
}
}
Vector3d baseline0 = Vector3d(0.,0.,0.);
BeamVis beammkn(baseline0, D_dish, thetaAntenne, phiAntenne, lambda0);
if (fgrectlobe) beammkn.setRectangularLobe(Lx,Ly);
if (fggaussbeam) beammkn.setGaussianLobe();
beammkn.AutoNormalizeBeam(mHeal);
*/
/* Aug 2018 : we should take into account the pixel size in the beam, when using this normalisation
factor with point sources ... -> beamnormfac * PixelSize = beamnormfac * 4 Pi / NbPixels() */
/*
double pixsize=1.;
{
SphereHEALPix<int_2> sphnp(mHeal);
pixsize=4.*M_PI/(double)sphnp.NbPixels();
}
double beamnormfac=beammkn.getBeamNormFactor()*pixsize;
cout << " ---> beamnormfac="<<beamnormfac<<" (PixelSize="<<pixsize<<" steradian)"<<endl;
// tm.Split("Done Normalisation computation");
double delta_pointing = latitude+shift_dec;
cout << "3- Computing visibilities for N= "<<dtk.NRows()<<" source positions - delta_pointing="
<<delta_pointing<<endl;
TArray< complex<double> > visiarr((sa_size_t)velev.size(), (sa_size_t)baselines.size());
for(sa_size_t jb=0; jb<(sa_size_t)baselines.size(); jb++) {
Vector3d baseline = baselines[jb];
if (prtlevel > 1)
cout << " Baseline["<<jb<<"] in tangent plane: " << baseline << endl;
BeamVis beammaker(baseline, D_dish, thetaAntenne, phiAntenne,lambda0);
if (fgrectlobe) {
beammaker.setRectangularLobe(Lx,Ly);
}
if (fggaussbeam) {
beammaker.setGaussianLobe();
}
beammaker.setBeamNormFactor(beamnormfac);
for(sa_size_t i=0; i<(sa_size_t)velev.size(); i++) {
double thetasrc=Angle(90.-velev[i],Angle::Degree).ToRadian();
double phisrcdeg=90.-vazim[i];
if (phisrcdeg<0.) phisrcdeg+=360.;
double phisrc=Angle(phisrcdeg,Angle::Degree).ToRadian();
visiarr(i, jb) = beammaker.Value(thetasrc, phisrc);
}
} // end of loop over baselines
cout << "4- Filling output DataTable ... "<<endl;
DataTable dtvis(128);
dtvis.AddDateTimeColumn("datetime");
dtvis.AddDoubleColumn("timesec");
dtvis.AddDoubleColumn("src_azimuth");
dtvis.AddDoubleColumn("src_elev");
for(size_t jb=0; jb<baselines.size(); jb++) {
char vcname[64];
sprintf(vcname,"V%d_%d",baselines_pair[jb].first,baselines_pair[jb].second);
dtvis.AddDoubleComplexColumn(vcname);
}
dtvis.Info().Comment()="Visibility from Satellite track DataTable";
DataTableRow row = dtvis.EmptyRow();
for(size_t i=0; i<velev.size(); i++) {
row[0] = vdate[i]; row[1] = vtms[i];
row[2] = vazim[i]; row[3] = velev[i];
for(size_t jb=0; jb<baselines.size(); jb++)
row[4+jb] = visiarr((sa_size_t)i, (sa_size_t)jb);
dtvis.AddRow(row);
}
dtvis.SetShowMinMaxFlag(true);
cout<<dtvis;
if (!fgoutfits) { // Sauvegarde du tableau de visibilites au format PPF
cout << "5- ----- Saving visibility DataTable to PPF file "<<outfilename<<endl;
POutPersist pos(outfilename);
pos << dtvis;
}
else { // // Sauvegarde du tableau de visibilites au format FITS
cout << "5- ----- Saving visibility DataTable FITS file "<<outfilename<<endl;
FitsInOutFile fos(outfilename,FitsInOutFile::Fits_Create);
fos << dtvis;
}
*/
} // End of try bloc
catch (PThrowable & exc) {
cerr << " trk_ac_fit.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 "
<< " - what()= " << e.what() << endl;
rc = 98;
}
catch (...) {
; cerr << " trk_ac_fit.cc: some other exception (...) was caught ! " << endl;
rc = 97;
}
cout << " ------------ END OF trk_ac_fit.cc (Rc= " << rc << ") ------------ " << endl;
return rc;
}
//-----------------------------------------------------------------------------------------
//-----------------------------------------------------------------------------------------
/* --Fonction-- */
/*
void setup_paon4()
{
cout << " --- Setting up PAON4 latitude, dish positions, dish size and frequency band ..." << endl;
dish_pos.clear();
Dgeom=paon4_Dgeom;
effD=paon4_Deff;
D_dish=paon4_D_dish;
for (int i=0; i<PAON4NBANT; i++)
dish_pos.push_back(Vector3d(paon4_posx[i], paon4_posy[i], paon4_posz[i]));
latitude=paon4_latitude;
freq0=1400.;
double clight = PhysQty::c().SIValue();
lambda0 = clight/(freq0*1.e6);
return;
}
*/
/* --Fonction-- */
/*
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"
<< " 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;
cout << " -paon4 : setup array as PAON4 \n"
<< " -lat latitude : define array latitude (in degree) \n"
<< " -D dish_diameter : define effective dish diameter (in meter D*eff) \n"
<< " -freq f0_MHz : define observation frequency , in MHz \n"
<< " -lambda lambda0_meter : define observation wavelength , in meters \n"
<< " -sdec angle : shift in antenna pointing in declination, in decimal degrees \n"
<< " as offset with respect to zenith, in the N-S plane, positive toward north \n"
<< " -in InputTrackPPFFile : Input source track file name (PPF datatable) \n"
<< " -out OutVisi_FileName : Output file (PPF or FITS) for the visibility DataTable file (def trkvis.ppf) \n"
<< " Specify path/name.fits ( filetype = fits) to write output array as a fits file \n"
<< " -prt PrintLevel: specify print level \n" << endl;
return 1;
}
infilename = "";
outfilename = "trkvis.ppf";
fgoutfits = false;
prtlevel=0;
double clight = PhysQty::c().SIValue();
freq0 = 1400.;
lambda0 = clight/(freq0*1.e6);
vector<string> lastargs;
while (narg>1) {
string fbo = arg[1];
if (fbo=="-paon4") { // Paon4 setup
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; }
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; }
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; }
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; }
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; }
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; }
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; }
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; }
prtlevel=atoi(arg[2]); arg+=2; narg-=2; lastargs.clear();
}
else { arg++; narg--; lastargs.push_back(fbo); }
}
if (lastargs.size() > 1) {
dish_pos.clear();
for (size_t i=0; i<lastargs.size(); i++) {
double x,y,z;
sscanf(lastargs[i].c_str(),"%lg,%lg,%lg",&x,&y,&z);
dish_pos.push_back(Vector3d(x,y,z));
}
}
if (dish_pos.size() < 2) {
cout << "trk_ac_fit/decode_args()/ERROR specify at least two dish positions ... " << endl;
return 2;
}
cout << " ---- trk_ac_fit/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++) {
cout << " Dish["<<i+1<<"] @ (x,y,z) meters (" << dish_pos[i].X()<<" , "<< dish_pos[i].Y()<<" , " << dish_pos[i].Z()<<")"<<endl;
}
//--- The auto-correlation baseline
baselines.push_back(dish_pos[0]-dish_pos[0]);
baselines_pair.push_back( pair<int, int>((int)1, (int)1) );
for (size_t i=0; i<dish_pos.size(); i++) {
for (size_t j=i+1; j<dish_pos.size(); j++) {
baselines.push_back(dish_pos[j]-dish_pos[i]);
baselines_pair.push_back( pair<int, int>((int)i+1, (int)j+1) );
}
}
cout << " --- Number of Baselines="<<baselines.size()<<" (X +toward East, Y +toward North, Z +toward Zenith)" << endl;
for (size_t k=0; k<baselines.size(); k++) {
cout << " Baseline["<<k+1<<"]->("<<baselines_pair[k].first<<","<<baselines_pair[k].second<<")= (dx,dy,dz)= ("
<< baselines[k].X()<<" , "<< baselines[k].Y()<<" , " << baselines[k].Z()<<") meters"<<endl;
}
size_t ll = outfilename.length();
size_t pp = outfilename.rfind('.');
if ((pp<ll)&&(outfilename.substr(pp)==".fits")) fgoutfits=true;
else fgoutfits=false;
cout << " Output file name: "<< outfilename << (fgoutfits?" [ FITS format ]":" [ PPF format]") << endl;
return 0;
}
*/
......@@ -6,8 +6,10 @@
avant l'appel a anasat.cc
FREQ_MIN , FREQ_MAX : definition de la bande de frequence en MHz , 1275-1282 pour Galileo
TIM_START , TIM_END : bornes en temps (en minutes) pour les calculs de normalisations
NORMFAC1,2,3,4 : argument optionnel - utiliser comme normalisation des 4 antennes voies 1H
c++args $FREQ_MIN $FREQ_MAX $TIM_START $TIM_END
c++execfrf anasat.cc
*/
if (args.size() < 4) {
cout << " anasat.cc / Missing arguments FREQ_MIN , FREQ_MAX , TIM_START , TIM_END"<<endl
......@@ -21,6 +23,16 @@ string FREQ_MAX = args[1];
string TIM_START = args[2];
string TIM_END = args[3];
bool fgnormarg=false;
double normarg[4]={1.,1.,1.,1.};
if (args.size() > 4) {
sscanf(args[4].c_str(), "%lg,%lg,%lg,%lg",normarg,normarg+1,normarg+2,normarg+3);
cout << " anasat.cc : Decoded normalisation argument: ";
for(size_t ii=0; ii<4; ii++) cout<<normarg[ii]<<" , ";
cout << endl;
fgnormarg=true;
}
// Frequency Binning
// int JMIN=51, JMAX=65;
// int JMIN=102, JMAX=132;
......@@ -63,23 +75,37 @@ for(sa_size_t i=0; i<V11.Size(); i++) {
V34(i)+=TFM_3H4H(i,j);
}
}
double fnorm=1./(double)(JMAX-JMIN+1);
complex<double> znorm=complex<double>(fnorm,0.);
V11 /= fnorm; V22 /= fnorm;
V33 /= fnorm; V44 /= fnorm;
V12 /= znorm; V13 /= fnorm; V14 /= znorm;
V23 /= znorm; V24 /= fnorm; V34 /= znorm;
V11 *= fnorm; V22 *= fnorm;
V33 *= fnorm; V44 *= fnorm;
V12 *= znorm; V13 *= fnorm; V14 *= znorm;
V23 *= znorm; V24 *= fnorm; V34 *= znorm;
double min,max;
vector<double> vfmax;
V11(Range(ITMA,ITMB)).MinMax(min,max);
cout << " V11->Max="<<max<<endl;
V22(Range(ITMA,ITMB)).MinMax(min,max); vfmax.push_back(max);
cout << " V22->Max="<<max<<endl;
V33(Range(ITMA,ITMB)).MinMax(min,max); vfmax.push_back(max);
cout << " V33->Max="<<max<<endl;
V44(Range(ITMA,ITMB)).MinMax(min,max); vfmax.push_back(max);
cout << " V44->Max="<<max<<endl;
vector<double> vfnorm;
V11(Range(ITMA,ITMB)).MinMax(min,max); vfnorm.push_back(max);
cout << " V11->Max="<<max<<endl; V11 /= max;
V22(Range(ITMA,ITMB)).MinMax(min,max); vfnorm.push_back(max);
cout << " V22->Max="<<max<<endl; V22 /= max;
V33(Range(ITMA,ITMB)).MinMax(min,max); vfnorm.push_back(max);
cout << " V33->Max="<<max<<endl; V33 /= max;
V44(Range(ITMA,ITMB)).MinMax(min,max); vfnorm.push_back(max);
cout << " V44->Max="<<max<<endl; V44 /= max;
if (fgnormarg)
for(size_t ii=0; ii<4; ii++) vfnorm.push_back(normarg[ii]);
else
for(size_t ii=0; ii<4; ii++) vfnorm.push_back(vfmax[ii]);
V11 /= vfnorm[0]; V22 /= vfnorm[1];
V33 /= vfnorm[2]; V44 /= vfnorm[3];
V12 /= complex<double>(sqrt(vfnorm[0]*vfnorm[1]));
V13 /= complex<double>(sqrt(vfnorm[0]*vfnorm[2]));
V14 /= complex<double>(sqrt(vfnorm[0]*vfnorm[3]));
......
......@@ -45,3 +45,72 @@
~/Work/Jiao/JSkyMap/Objs/trk2vis -paon4 -in trk_B01_41550.ppf -freq 1285 -sdec 9 -out trkvis_B01_41550.ppf
~/Work/Jiao/JSkyMap/Objs/trk2vis -paon4 -in trk_B02_40890.ppf -freq 1285 -sdec 9 -out trkvis_B02_40890.ppf
~/Work/Jiao/JSkyMap/Objs/trk2vis -paon4 -in trk_B04_40889.ppf -freq 1285 -sdec 11.9 -out trkvis_B04_40889.ppf
########################################
### Analyse de A21_14oct18 B04_26oct18 B12_6dec18
### A21_14oct18
~/Work/GIT_BAORadio/AnaPAON4/Satellites/Objs/predictsatsgp4 -T "2018/10/14 23:00:00" -H 0.,78.58 -K 43057 TLE_20181123/*
~/Work/GIT_BAORadio/AnaPAON4/Satellites/Objs/predictsatsgp4 -T "2018/10/15 1:33:00" -H 0.,78.58 -K 43055 TLE_20181123/*
~/Work/GIT_BAORadio/AnaPAON4/Satellites/Objs/trk2dt trk_43057_20181014.txt trk_A21_43057.ppf 2018/10/14
~/Work/GIT_BAORadio/AnaPAON4/Satellites/Objs/trk2dt trk_43055_20181015.txt trk_A21_43055.ppf 2018/10/14