Skip to content

GitLab

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

Ajout script pour analyse satellites PAON4, Reza 4/12/2018

parent a563df01
Hide whitespace changes
Inline Side-by-side
Showing with 365 additions and 0 deletions
Scripts/anasat.cc 0 → 100644
View file @ bafca665
/*
Analyse PAON4 , Premier essais sur les signaux des satellites
R. Ansari, Novembre 2018
Code destine a etre utilise avec anasat.pic
Il faut definir les 4 arguments suivants au niveau de l'interpreteur par c++args
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
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
<< " Define FREQ_MIN , FREQ_MAX , TIM_START , TIM_END "<<endl
<< " Use c++args $FREQ_MIN $FREQ_MAX $TIM_START $TIM_END"<<endl
<< " before calling anasat.cc through c++exec " << endl;
return 9;
}
string FREQ_MIN = args[0];
string FREQ_MAX = args[1];
string TIM_START = args[2];
string TIM_END = args[3];
// Frequency Binning
// int JMIN=51, JMAX=65;
// int JMIN=102, JMAX=132;
double delta_freq=250./(double)TFM_1H.SizeY();
sa_size_t JMIN = (atof(FREQ_MIN.c_str())-1250.+0.5)/delta_freq;
sa_size_t JMAX = (atof(FREQ_MAX.c_str())-1250.+0.5)/delta_freq;
cout << "anasat.cc/Info: FREQ_MIN,MAX="<<FREQ_MIN<<","<<FREQ_MAX<<" -> JMIN,MAX="<<JMIN<<","<<JMAX<<endl;
// Normalisation Time-range
// sa_size_t ITMA=180, ITMB=220;
// sa_size_t ITMA=360, ITMB=440; /* B01 41859 */
// sa_size_t ITMA=1010, ITMB=1070; /* B04 40889 */
double tm_start=TimeVec(0);
double tm_end=TimeVec(TimeVec.Size()-1);
double delta_temps=(tm_end-tm_start)/(double)TimeVec.Size();
sa_size_t ITMA=(atof(TIM_START.c_str())*60.-tm_start+0.5)/delta_temps;
sa_size_t ITMB=(atof(TIM_END.c_str())*60.-tm_start+0.5)/delta_temps;
cout << "anasat.cc/Info: TIM_START,END="<<TIM_START<<","<<TIM_END<<" -> ITMA,B="<<ITMA<<","<<ITMB<<endl;
sa_size_t SZX=TFM_1H.SizeX();
Vector V11(SZX); V11=0.; Vector V22(SZX); V22=0.;
Vector V33(SZX); V33=0.; Vector V44(SZX); V44=0.;
TVector< complex<double> > V12(SZX); V12= complex<double>(0.,0.);
TVector< complex<double> > V13(SZX); V13= complex<double>(0.,0.);
TVector< complex<double> > V14(SZX); V14= complex<double>(0.,0.);
TVector< complex<double> > V23(SZX); V23= complex<double>(0.,0.);
TVector< complex<double> > V24(SZX); V24= complex<double>(0.,0.);
TVector< complex<double> > V34(SZX); V34= complex<double>(0.,0.);
for(sa_size_t i=0; i<V11.Size(); i++) {
for(sa_size_t j=JMIN; j<=JMAX; j++) {
V11(i)+=TFM_1H(i,j); V22(i)+=TFM_2H(i,j);
V33(i)+=TFM_3H(i,j); V44(i)+=TFM_4H(i,j);
V12(i)+=TFM_1H2H(i,j);
V13(i)+=TFM_1H3H(i,j);
V14(i)+=TFM_1H4H(i,j);
V23(i)+=TFM_2H3H(i,j);
V24(i)+=TFM_2H4H(i,j);
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;
double min,max;
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;
V12 /= complex<double>(sqrt(vfnorm[0]*vfnorm[1]));
V13 /= complex<double>(sqrt(vfnorm[0]*vfnorm[2]));
V14 /= complex<double>(sqrt(vfnorm[0]*vfnorm[3]));
V23 /= complex<double>(sqrt(vfnorm[1]*vfnorm[2]));
V24 /= complex<double>(sqrt(vfnorm[1]*vfnorm[3]));
V34 /= complex<double>(sqrt(vfnorm[2]*vfnorm[3]));
cout << " Finished normalisation --- Creating DataTable ..."<<endl;
DataTable dtvis(128);
dtvis.AddDoubleColumn("timesec");
dtvis.AddDoubleColumn("ra_hour");
dtvis.AddDoubleColumn("V11");
dtvis.AddDoubleColumn("V22");
dtvis.AddDoubleColumn("V33");
dtvis.AddDoubleColumn("V44");
dtvis.AddDoubleComplexColumn("V12");
dtvis.AddDoubleComplexColumn("V13");
dtvis.AddDoubleComplexColumn("V14");
dtvis.AddDoubleComplexColumn("V23");
dtvis.AddDoubleComplexColumn("V24");
dtvis.AddDoubleComplexColumn("V34");
DataTableRow row = dtvis.EmptyRow();
for(size_t i=0; i<(size_t)SZX; i++) {
row[0] = TimeVec(i);
row[1] = RAVec(i);
row[2] = V11(i); row[3] = V22(i);
row[4] = V33(i); row[5] = V44(i);
row[6] = V12(i); row[7] = V13(i);
row[8] = V14(i); row[9] = V23(i);
row[10] = V24(i); row[11] = V34(i);
dtvis.AddRow(row);
}
dtvis.SetShowMinMaxFlag(true);
cout<<dtvis;
KeepObj(dtvis);
Scripts/anasat.pic 0 → 100644
View file @ bafca665
##########################################################################
####### Analyse PAON4 , Premier essais sur les signaux des satellites
####### R. Ansari, Novembre 2018
## Analyse / affichage a partir des cartes TFM Auto et Cross-Correlation
## et DataTable visibilites fait par le programme JSkyMap/trk2vis.cc
## Ce script utilise anasat.cc
##########################################################################
setaxesatt 'font=helvetica,bold,16 fixedfontsize minorticks'
###########################################
###### Part 1 : Analyse B01
###########################################
delobjs *
openppf tfm_B01_20oct18.ppf
FREQ_MIN = 1275
FREQ_MAX = 1282
TIM_START = 900
TIM_END = 950
c++args $FREQ_MIN $FREQ_MAX $TIM_START $TIM_END
c++execfrf anasat.cc
openppf trk_B01_41859.ppf
rename trk_B01_41859 trk
openppf trkvis_B01_41859.ppf
rename trkvis_B01_41859 trkvis
openppf trk_B01_41550.ppf
rename trk_B01_41550 trk2
openppf trkvis_B01_41550.ppf
rename trkvis_B01_41550 trkvis2
plot2d trk (azimuth<100?azimuth+360:azimuth) 90-elevation 1 'line=solid,2 nsta cpts red'
plot2d trk2 (azimuth<100?azimuth+360:azimuth) 90-elevation 1 'line=solid,2 same nsta cpts magenta'
addline 260 11 420 11 'line=solid,4 black'
# addline 260 7 420 7 'line=solid,4 black'
settitle 'B01 Galileo 41859 (red) 41550 (magenta) ' ' ' 'font=helvetica,bold,18 black'
setaxelabels 'Azimuth (degree)' 'Zenith Angle (=90-Elevation) (degree)' 'font=helvetica,bolditalic,18 black'
# addcirc 360. 9 3 'line=solid,4 black'
# addoval 2.19 56.3 7 4 'black line=solid,4'
# addoval 2.19 58.5 7 4 'black line=solid,4'
plot2d trk (azimuth<100?azimuth+360:azimuth) 90-elevation 1 'line=solid,2 nsta cpts red'
addcirc 360. 12 3 'line=solid,4 black'
plot2d dtvis timesec/60 V11 1 'cpts black nsta notit'
plot2d dtvis timesec/60 V22 1 'cpts navyblue same nsta'
plot2d dtvis timesec/60 V33 1 'cpts blue same nsta'
plot2d dtvis timesec/60 V44 1 'cpts cyan same nsta'
plot2d trkvis timesec/60 V1_1.real*1.1e3 1 'same red line=dashed,2 cpts nsta'
plot2d trkvis2 timesec/60 V1_1.real*1.5e3 1 'same magenta line=dashed,2 cpts nsta'
plot2d dtvis timesec/60 V12.real 1 'cpts line=solid,2 black nsta notit'
plot2d dtvis timesec/60 V12.imag 1 'cpts same blue line=solid,2 marker=box,5 nsta notit'
settitle '1H-2H real/imag , dashed=predited' ' ' 'font=helvetica,bolditalic,16 black'
# plot2d trkvis timesec/60 V1_2.real*1.1e3 1 'same red line=dashed,2 cpts nsta'
plot2d trkvis timesec/60 -V1_2.imag*1.1e3 1 'same orange line=dashed,2 cpts nsta'
plot2d trkvis2 timesec/60 V1_2.real*1.1e3 1 'same magenta line=dashed,2 cpts nsta'
plot2d trkvis2 timesec/60 -V1_2.imag*1.1e3 1 'same violetred line=dashed,2 cpts nsta'
phi = 0.65*Pi
cp = cos(phi)
sp = sin(phi)
A = 0.95e3
set RPHI (V1_2.real*$cp-V1_2.imag*$sp)*$A
set IPHI -(V1_2.imag*$cp+V1_2.real*$sp)*$A
plot2d dtvis timesec/60 V12.real 1 'xylimits=900,950,-1,1 cpts line=solid,2 marker=fbox,7 black nsta notit'
plot2d dtvis timesec/60 V12.imag 1 'cpts same blue line=solid,2 marker=fbox,7 nsta notit'
plot2d trkvis timesec/60 $RPHI 1 'same red line=dashed,2 marker=circle,5 cpts nsta'
plot2d trkvis timesec/60 $IPHI 1 'same magenta line=dashed,2 marker=circle,5 cpts nsta'
settitle '1H-2H real/imag , dashed=predited B01 Galileo 41859' ' ' 'font=helvetica,bolditalic,16 black'
plot2d trkvis2 timesec/60 $RPHI 1 'same red line=dashed,2 marker=circle,5 cpts nsta'
plot2d trkvis2 timesec/60 $IPHI 1 'same magenta line=dashed,2 marker=circle,5 cpts nsta'
settitle 'B01 Galileo 41859 1H2H CrossCorr Predicted Satellite' ' ' 'font=helvetica,bold,18 black'
setaxelabels 'Time (min)' 'CrossCor Real/Imag' 'font=helvetica,bolditalic,18 black'
phi = 0.12*Pi
cp = cos(phi)
sp = sin(phi)
A = 0.95e3
set RPHI (V3_4.real*$cp-V3_4.imag*$sp)*$A
set IPHI -(V3_4.imag*$cp+V3_4.real*$sp)*$A
plot2d dtvis timesec/60 V34.real 1 'xylimits=900,950,-1,1 cpts line=solid,2 marker=fbox,7 black nsta notit'
plot2d dtvis timesec/60 V34.imag 1 'cpts same blue line=solid,2 marker=fbox,7 nsta notit'
plot2d trkvis timesec/60 $RPHI 1 'same red line=dashed,2 marker=circle,5 cpts nsta'
plot2d trkvis timesec/60 $IPHI 1 'same magenta line=dashed,2 marker=circle,5 cpts nsta'
settitle 'B01 Galileo 41859 3H4H CrossCorr Predicted Satellite' ' ' 'font=helvetica,bold,18 black'
setaxelabels 'Time (min)' 'CrossCor Real/Imag' 'font=helvetica,bolditalic,18 black'
plot2d dtvis timesec/60 V23.real 1 'cpts line=solid,2 black nsta notit'
plot2d dtvis timesec/60 V23.imag 1 'cpts same blue line=solid,2 marker=box,5 nsta notit'
settitle '2H-3H real/imag , dashed=predited' ' ' 'font=helvetica,bolditalic,16 black'
plot2d trkvis timesec/60 V2_3.real*1.1e3 1 'same red line=dashed,2 cpts nsta'
plot2d trkvis timesec/60 -V2_3.imag*1.1e3 1 'same orange line=dashed,2 cpts nsta'
plot2d trkvis2 timesec/60 V2_3.real*1.1e3 1 'same magenta line=dashed,2 cpts nsta'
plot2d trkvis2 timesec/60 -V2_3.imag*1.1e3 1 'same violetred line=dashed,2 cpts nsta'
plot2d dtvis timesec/60 V34.real 1 'cpts line=solid,2 black nsta notit'
plot2d dtvis timesec/60 V34.imag 1 'cpts same blue line=solid,2 marker=box,5 nsta notit'
settitle '3H-4H real/imag , dashed=predited' ' ' 'font=helvetica,bolditalic,16 black'
plot2d trkvis timesec/60 V3_4.real*1.1e3 1 'same red line=dashed,2 cpts nsta'
plot2d trkvis timesec/60 -V3_4.imag*1.1e3 1 'same orange line=dashed,2 cpts nsta'
plot2d trkvis2 timesec/60 V3_4.real*1.1e3 1 'same magenta line=dashed,2 cpts nsta'
plot2d trkvis2 timesec/60 -V3_4.imag*1.1e3 1 'same violetred line=dashed,2 cpts nsta'
disp TFM_1H 'colbr128 stdaxes wbgcol=white'
###########################################
###### Part 1 : Analyse B04
###########################################
delobjs *
openppf tfm_B01_20oct18.ppf
FREQ_MIN = 1275
FREQ_MAX = 1282
TIM_START = 900
TIM_END = 950
c++args $FREQ_MIN $FREQ_MAX $TIM_START $TIM_END
c++execfrf anasat.cc
delobjs *
openppf tfm_B04_26oct18.ppf
FREQ_MIN = 1275
FREQ_MAX = 1282
TIM_START = 1420
TIM_END = 1500
c++args $FREQ_MIN $FREQ_MAX $TIM_START $TIM_END
c++execfrf anasat.cc
openppf trk_B04_40889.ppf
rename trk_B04_40889 trk
openppf trkvis_B04_40889.ppf
rename trkvis_B04_40889 trkvis
plot2d trk (azimuth<100?azimuth+360:azimuth) 90-elevation 1 'line=solid,2 nsta cpts red'
addline 300 14 420 14 'line=solid,4 black'
# addline 300 9 420 9 'line=solid,4 black'
settitle 'B04 Galileo 40889' ' ' 'font=helvetica,bold,18 black'
setaxelabels 'Azimuth (degree)' 'Elevation (degree)' 'font=helvetica,bolditalic,18 black'
plot2d dtvis timesec/60 V11 1 'cpts black nsta notit xylimits=1420,1500,0.,1.1'
plot2d dtvis timesec/60 V22 1 'cpts navyblue same nsta'
plot2d dtvis timesec/60 V33 1 'cpts blue same nsta'
plot2d dtvis timesec/60 V44 1 'cpts cyan same nsta'
plot2d trkvis timesec/60 V1_1.real*1.3e3+0.05 1 'same red line=dashed,2 cpts nsta'
settitle 'B04 Galileo 40889 4 HH AutoCorr + Predicted Satellite (red)' ' ' 'font=helvetica,bold,18 black'
setaxelabels 'Time (min)' 'Normalized AutoCor' 'font=helvetica,bolditalic,18 black'
phi = 0.65*Pi
cp = cos(phi)
sp = sin(phi)
A = 1.2e3
set RPHI (V1_2.real*$cp-V1_2.imag*$sp)*$A
set IPHI -(V1_2.imag*$cp+V1_2.real*$sp)*$A
plot2d dtvis timesec/60 V12.real 1 'xylimits=1430,1490,-1,1 cpts line=solid,2 marker=fbox,7 black nsta notit'
plot2d dtvis timesec/60 V12.imag 1 'cpts same blue line=solid,2 marker=fbox,7 nsta notit'
plot2d trkvis timesec/60 $RPHI 1 'same red line=dashed,2 marker=circle,5 cpts nsta'
plot2d trkvis timesec/60 $IPHI 1 'same magenta line=dashed,2 marker=circle,5 cpts nsta'
settitle 'B04 Galileo 40889 1H2H CrossCorr Predicted Satellite' ' ' 'font=helvetica,bold,18 black'
setaxelabels 'Time (min)' 'CrossCor Real/Imag' 'font=helvetica,bolditalic,18 black'
phi = 0.12*Pi
cp = cos(phi)
sp = sin(phi)
A = 1.2e3
set RPHI (V3_4.real*$cp-V3_4.imag*$sp)*$A
set IPHI -(V3_4.imag*$cp+V3_4.real*$sp)*$A
plot2d dtvis timesec/60 V34.real 1 'xylimits=1430,1490,-1,1 cpts line=solid,2 marker=fbox,7 black nsta notit'
plot2d dtvis timesec/60 V34.imag 1 'cpts same blue line=solid,2 marker=fbox,7 nsta notit'
plot2d trkvis timesec/60 $RPHI 1 'same red line=dashed,2 marker=circle,5 cpts nsta'
plot2d trkvis timesec/60 $IPHI 1 'same magenta line=dashed,2 marker=circle,5 cpts nsta'
settitle 'B04 Galileo 40889 3H4H CrossCorr Predicted Satellite' ' ' 'font=helvetica,bold,18 black'
setaxelabels 'Time (min)' 'CrossCor Real/Imag' 'font=helvetica,bolditalic,18 black'
objaoper TFM_2H yline 1048
plot2d yline_1048 1250+n*250/1024 val/26e6 1 'cpts magenta nsta notit line=solid,2 xylimits=1250,1325,0.,1.1'
settitle 'B04 Galileo 40889 Freq. Spectrum' ' ' 'font=helvetica,bold,18 black'
setaxelabels 'Freq (MHz)' 'Power (au)' 'font=helvetica,bolditalic,18 black'
Scripts/p4sat.cmd 0 → 100644
View file @ bafca665
######################################################################
####### Exemples de commandes pour analyses des satellites ds PAON4
####### Analyse PAON4 ,
######## 1ere version , R. Ansari, Novembre 2018
######################################################################
#####################
## A/ On utilise le programme predictsatsgp4 pour generer des fichiers de traces (position en fonction du temps)
## pour les satellites
#####################
###### B01 , 20 Oct 2018
# Satellites vers 15:21 UT de B01 Galileo 41859
~/Work/GIT_BAORadio/AnaPAON4/Satellites/Objs/predictsatsgp4 -T "2018/10/20 15:21:18" -H 0.,81.1 -K 41859 -S 10 -s 1,0,0,15 -p 0 TLE_20181123/galileo.txt
# Satellites vers 23:40 UT de B01 Galileo 41550
~/Work/GIT_BAORadio/AnaPAON4/Satellites/Objs/predictsatsgp4 -T "2018/10/20 23:39:59.0" -H 0.,81.1 -K 41550 -S 10 -s 1,0,0,15 -p 0 TLE_20181123/galileo.txt
###### B02 , 21 Oct 2018
# Satellites vers 22:40 UT de B02 Galileo 40890
~/Work/GIT_BAORadio/AnaPAON4/Satellites/Objs/predictsatsgp4 -T "2018/10/21 22:02:0.0" -H 0.,82.1 -K 40890 -S 10 -s 1,0,0,15 -p 0 TLE_20181123/galileo.txt
###### B04 , 26 Oct 2018
# Satellites vers 27/10/2018 00:18:00 UT de B04
~/Work/GIT_BAORadio/AnaPAON4/Satellites/Objs/predictsatsgp4 -T "2018/10/27 00:18:0.0" -H 0.,78.1 -K 40889 -s 1,0,0,15 -p 0 galileo.txt
#####################
## B/ On convertit les fichiers de trace au format texte en DataTable PPF ou FITS (bin-table) avec trk2dt
#####################
#### Converting to PPF DataTable track files
~/Work/GIT_BAORadio/AnaPAON4/Satellites/Objs/trk2dt trk_41859_20181020.txt trk_B01_41859.ppf
~/Work/GIT_BAORadio/AnaPAON4/Satellites/Objs/trk2dt trk_41550_20181020.txt trk_B01_41550.ppf
~/Work/GIT_BAORadio/AnaPAON4/Satellites/Objs/trk2dt trk_40890_20181021.txt trk_B02_40890.ppf
~/Work/GIT_BAORadio/AnaPAON4/Satellites/Objs/trk2dt trk_40889_20181027.txt trk_B04_40889.ppf
#####################
## C/ On calcule les visibilites PAON4 a partir des fichiers de traces format PPF
## a l'aide du programme JSkyMap trk2vis
#####################
#### Computing visibility DataTable from track files
~/Work/Jiao/JSkyMap/Objs/trk2vis -paon4 -in trk_B01_41859.ppf -freq 1285 -sdec 9 -out trkvis_B01_41859.ppf
~/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
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