Commit 1b8f91a6 authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

Modif classe CxBeam et ajout classes My6CxSignalsB / My6CxGenXi2B pour...

Modif classe CxBeam et ajout classes My6CxSignalsB / My6CxGenXi2B pour ajustement des baselines, Reza 13/02/2019
parent 12e96e08
......@@ -26,6 +26,12 @@ using namespace SOPHYA;
//-------------------------------------------------------------------------
class ACBeam { // Auto-correlation beam
public:
// Default constructor , D=5m diameter beam lambda=0.21
ACBeam()
: fggauss_(false), D_(5.), lambda_(0.21), DoL_(5./0.21),
bdir_(LongitudeLatitude(0., 0.)), bdir_theta_(0.), bdir_phi_(0.), NormFac_(1.)
{
}
// constructor: define a lobe to a circular dish of diameter D pointing to the direction (theta_b,phi_b)
ACBeam(double D, double theta_b=0., double phi_b=0., double lambda=0.21)
: fggauss_(false), D_(D), lambda_(lambda), DoL_(D/lambda),
......
......@@ -14,11 +14,20 @@ using namespace SOPHYA;
class CxBeam { // Auto-correlation beam
public:
// Default constructor : an auto-cor beam for default ACBeam
CxBeam()
: a1_(ACBeam()), a2_(ACBeam()), base_baseline_(Vector3d(0.,0.,0.)), baseline_(Vector3d(0.,0.,0.))
{
if (fabs(a1_.getLambda()-a2_.getLambda())>1e-9)
throw ParmError("CxBeam::CxBeam() - not same lambda in the two different beams !");
facphase_=2.*M_PI/a1_.getLambda();
}
CxBeam(ACBeam & a1, ACBeam & a2, Vector3d & baseline)
: a1_(a1), a2_(a2), base_baseline_(baseline), baseline_(baseline)
{
if (fabs(a1.getLambda()-a2.getLambda())>1e-9)
throw ParmError("CxBeam::CxBeam() - not same lambda in the two different beams !");
throw ParmError("CxBeam::CxBeam(a1,a2,baseline) - not same lambda in the two different beams !");
facphase_=2.*M_PI/a1.getLambda();
}
CxBeam(CxBeam const & c)
......
/* PAON4 analysis software
Class for fitting baselines and instrumental phases using the 6 baselines
R. Ansari, February 2019 */
#ifndef GCXFITBASELINE_H_SEEN
#define GCXFITBASELINE_H_SEEN
#include "generalfit.h"
#include "simplex.h"
#include "sunitpcst.h"
#include "cxbeam.h"
class My6CxSignalsB {
public:
My6CxSignalsB(vector< vector<double> > & v_time_data_, vector< vector< vector< complex<double> > > > & vv_cxdata_,
vector< vector< vector<double> > > & vv_cxerr_, // vector<double> & v_freqs_,
vector< SLinInterp1D > & v_interp_elev_, vector< SLinInterp1D > & v_interp_sazim_,
vector< vector<CxBeam> > & v_beams_, vector< vector<double> > & v_A_) // vector< vector< complex<double> > > & v_B_
: n_cxcor(6), v_time_data(v_time_data_), vv_cxdata(vv_cxdata_), vv_cxerr(vv_cxerr_), // v_freqs(v_freqs_),
v_interp_elev(v_interp_elev_), v_interp_sazim(v_interp_sazim_),
v_beams(v_beams_), v_A(v_A_)
{
//DBG cout << " My6CxSignalsB() "<< endl;
}
// get the j-th signal for cross-correlation kcx (0<=kcx<=5)
virtual int getDataSignal(size_t j, size_t kcx, vector< complex<double> > & sig)
{
//DEL sig.resize(v_time_data[j].size());
vector< complex<double> > & vcxdata = vv_cxdata[j][kcx];
sig = vcxdata;
//DEL for(size_t k=0; k<sig.size(); k++) sig[k]=vcxdata[k];
return (int)sig.size();
}
// get the j-th signal for cross-correlation kcx (0<=kcx<=5)
virtual TVector< complex<double> > getDataSignal(size_t j, size_t kcx)
{
vector< complex<double> >sig;
getDataSignal(j, kcx, sig);
TVector< complex<double> > rvs(sig);
//DEL for(size_t k=0; k<sig.size(); k++) rvs(k)=sig[k];
return rvs;
}
// get the beam for the j-th signal and the cross-correlation kcx (0<=kcx<=5)
inline CxBeam & getBeam(size_t j, size_t kcx)
{
vector<CxBeam> & beams = v_beams[j];
return beams[kcx];
}
// get the beam for the j-th signal and the cross-correlation kcx (0<=kcx<=5)
virtual int getExpectedSignal(size_t j, size_t kcx, vector< complex<double> > & sig, double phase, Vector3d & baseline_shift)
{
//DBG cout<<"*DBG*A*My6CxSignalsB::getExpectedSignal() j="<<j<<" size="<<v_time_data[j].size()<<" kcx="<<kcx<<endl;
sig.resize(v_time_data[j].size());
vector<double> & vtime = v_time_data[j];
// vector< complex<double> > & vcxdata = vv_cxdata[j][nac_];
// vector<double> & vcxerr = vv_cxerr[j][nac_];
vector<double> & lesA = v_A[j];
double A = lesA[kcx];
CxBeam beam = getBeam(j, kcx);
beam.ShiftBaseline(baseline_shift);
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 sazim=v_interp_sazim[j].YInterp(vtime[k]);
while (sazim > 360.) sazim -= 360.;
double phisrcdeg = 90.-sazim;
if (phisrcdeg < 0.) phisrcdeg+=360.;
complex<double> cxb=beam.Value(thetasrc, Angle(phisrcdeg,Angle::Degree).ToRadian());
// on fait conjugue complexe - car les donnees PAON4 sont Vi conj(Vj) , ici, on calcule conj(Vi) * Vj
sig[k] = complex<double>(A*cos(phase),A*sin(phase))*conj(cxb);
//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();
}
/* definition du tableau des parametres de fit
parm[0] = Phase_2
parm[1] = Phase_3
parm[2] = Phase_4
parm[3] = X_shift_2
parm[4] = Y_shift_2
parm[5] = Z_shift_2
parm[6] = X_shift_3
parm[7] = Y_shift_3
parm[8] = Z_shift_3
parm[9] = X_shift_4
parm[10] = Y_shift_4
parm[11] = Z_shift_4
*/
// determine la phase et le baseline_shift a partir du tableau parm, pour la ligne de base kcx (0<=kcx<=5)
void getFromFitParam(size_t kcx, const double* parm, double & phase, Vector3d & baseline_shift)
{
size_t na[6]={0,0,0,1,1,2};
size_t nb[6]={1,2,3,2,3,3};
size_t na0[6]={0,0,0,0,0,1};
size_t nb0[6]={0,1,2,1,2,2};
if (kcx<3) {
phase=parm[kcx];
baseline_shift=Vector3d(parm[3*kcx+3], parm[3*kcx+4], parm[3*kcx+5]);
}
else {
size_t ja=na0[kcx] , jb=nb0[kcx];
phase=parm[jb]-parm[ja];
baseline_shift=Vector3d(parm[3*jb+3]-parm[3*ja+3], parm[3*jb+4]-parm[3*ja+4], parm[3*jb+5]-parm[3*ja+5]);
}
return;
}
// get the beam for the j-th signal and the cross-correlation kcx (0<=kcx<=5)
virtual int getExpectedSignal(size_t j, size_t kcx, vector< complex<double> > & sig, const double* parm)
{
double phase;
Vector3d baseline_shift;
getFromFitParam(kcx, parm, phase, baseline_shift);
double A = parm[1+3*j];
complex<double> B = complex<double>(parm[2+3*j], parm[3+3*j]);
return getExpectedSignal(j, kcx, sig, phase, baseline_shift);
}
virtual TVector< complex<double> > getExpectedSignal(size_t j, size_t kcx, double phase, Vector3d & baseline_shift)
{
vector< complex<double> > sig;
//DBG cout << "*DBG*getExpectedSignal() : calling getExpectedSignal(j="<<j<<" ,sig,...)"<<endl;
//DBG debug_gacfit=1;
getExpectedSignal(j, kcx, sig, phase, baseline_shift);
TVector< complex<double> > rvs(sig);
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++) { // loop over trk data sets
for(size_t kcx=0; kcx<n_cxcor; kcx++) { // loop over the six cross-correlations
vector< complex<double> > sig;
npts += getExpectedSignal(j, kcx, sig, parm);
vector< complex<double> > & vcxdata = vv_cxdata[j][kcx];
vector<double> & verr = vv_cxerr[j][kcx];
for(size_t i=0; i<vcxdata.size(); i++) {
double xx=(sig[i].real()-vcxdata[i].real())/verr[i];
xi2 += (xx*xx);
xx=(sig[i].imag()-vcxdata[i].imag())/verr[i];
xi2 += (xx*xx);
}
}
}
npts *= 2;
return xi2;
}
size_t n_cxcor; // Number of Cx correlations ( 6 for PAON4 )
vector< vector<double> > & v_time_data;
vector< vector< vector< complex<double> > > > & vv_cxdata;
vector< vector< vector<double> > > & vv_cxerr;
vector< SLinInterp1D > & v_interp_elev;
vector< SLinInterp1D > & v_interp_sazim;
vector< vector<CxBeam> > & v_beams;
vector< vector<double> > & v_A; // vecteur des ampltudes pour chaque cross-correlation - determine par le fit d'avant
// vector< vector< complex<double> > > & v_B; // vecteur des offsets pour chaque cross-correlation - determine par le fit d'avant
};
//-----------------------------------------------------------------
class My6CxGenXi2B : public GeneralXi2, public My6CxSignalsB {
public:
My6CxGenXi2B(vector< vector<double> > & v_time_data_, vector< vector< vector< complex<double> > > > & vv_cxdata_,
vector< vector< vector<double> > > & vv_cxerr_, // vector<double> & v_freqs_,
vector< SLinInterp1D > & v_interp_elev_, vector< SLinInterp1D > & v_interp_sazim_,
vector< vector<CxBeam> > & v_beams_, vector< vector<double> > & v_A_)
: GeneralXi2(12) , // 12 = 3 phases + 3*3 position shifts
My6CxSignalsB(v_time_data_, vv_cxdata_, vv_cxerr_, v_interp_elev_, v_interp_sazim_, v_beams_, v_A_)
{
}
/* definition du tableau des parametres de fit
parm[0] = Phase_2 Phase Antenne 2 = DeltaPhi_12
parm[1] = Phase_3 Phase Antenne 3 = DeltaPhi_13
parm[2] = Phase_4 Phase Antenne 4 = DeltaPhi_14
parm[3] = X_shift_2 Shift (/nominal) X Antenne 2 ( DeltaX_12 )
parm[4] = Y_shift_2 Shift (/nominal) Y Antenne 2 ( DeltaY_12 )
parm[5] = Z_shift_2 Shift (/nominal) Z Antenne 2 ( DeltaZ_12 )
parm[6] = X_shift_3 Shift (/nominal) X Antenne 3 ( DeltaX_13 )
parm[7] = Y_shift_3 Shift (/nominal) Y Antenne 3 ( DeltaY_13 )
parm[8] = Z_shift_3 Shift (/nominal) Z Antenne 2 ( DeltaZ_13 )
parm[9] = X_shift_4 Shift (/nominal) X Antenne 4 ( DeltaX_14 )
parm[10] = Y_shift_4 Shift (/nominal) Y Antenne 4 ( DeltaY_14 )
parm[11] = Z_shift_4 Shift (/nominal) Z Antenne 4 ( DeltaZ_14 )
*/
virtual double Value(GeneralFitData& data, double* parm, int& ndataused)
{
return getXi2(parm, ndataused);
}
};
#endif
......@@ -38,7 +38,7 @@
#include "acbeam.h"
#include "gacfit.h"
#include "gcxfit.h"
#include "gcxfitbaseline.h"
using namespace std;
using namespace SOPHYA;
......
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