add code to use Jean-Eric's polynomial phase model + small changes to visiavg

p4phcor.cc

@@ -60,19 +60,29 @@ P4PhaseCor::P4PhaseCor(string const & filename)
   readFitParamFile(filename);
 }
 
+double P4PhaseCor::jec_phase_poly(double freq, double * p){
+  double  ret = 0.;
+  double  del = 1.;
+  for (int k=0 ; k<int(p[0])+1 ; k++,del *= freq-p[1] ){
+    ret += p[k+2]*del ;
+  }
+  return(ret);
+}
+
 void P4PhaseCor::readFitParamFile(string const & filename)
 {
   cout << "--- P4PhaseCor::readFitParamFile( filename= "<<filename<<" )"<<endl;
   std::ifstream is(filename.c_str(), std::ifstream::in);
+  
   string line;
   char mot[32], ftype[32], key[4], feed[3];
-  double p[3];
+  double p[12] ;
   while (!is.eof()) {
     line = "";
     getline(is, line);
     if ((is.good() || is.eof()) && (line.length()>0) && (line[0]!='#')) {
-      p[0]=p[1]=p[2]= -9.e-99;
-      sscanf(line.c_str(),"%s %s %lg %lg %lg", mot, ftype, p, p+1, p+2);
+      for (int k=0 ; k<12 ; k++) p[k] = -9e-99;
+      sscanf(line.c_str(),"%s %s %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg %lg", mot, ftype, p, p+1, p+2,p+3,p+4,p+5,p+6,p+7,p+8,p+9,p+10,p+11,p+12);
       key[0]=mot[0];  key[1]=mot[1];  key[2]=mot[2];  key[3]='\0';
       feed[0]=mot[3];  feed[1]=mot[4];  feed[2]='\0';
       if (strcmp(key,"Phi")!=0)  continue;
@@ -82,7 +92,7 @@ void P4PhaseCor::readFitParamFile(string const & filename)
 	cout<<"P4PhaseCor::readFitParamFile()/ERROR  bad feed name= "<<feed<<" -> fid="<<fid<<endl;
 	continue;
       }
-      if ((sftyp!="Sawtooth")&&(sftyp!="Linear"))  {
+      if ((sftyp!="Sawtooth")&&(sftyp!="Linear")&&(sftyp!="Polynomial"))  {
 	cout<<"P4PhaseCor::readFitParamFile()/ERROR  bad fitted function type= "<<sftyp<<endl;
 	continue;
       }
@@ -90,7 +100,8 @@ void P4PhaseCor::readFitParamFile(string const & filename)
       for (sa_size_t k=0; k<phases.Size(); k++) {
 	double freq = P4FreqBand::getP4Frequency((int)k);
 	if (sftyp=="Linear")  phases(k) = jec_phase_linear(freq, p);
-	else phases(k) = jec_phase_sawtooth(freq, p);
+	else if (sftyp=="Sawtooth") phases(k) = jec_phase_sawtooth(freq, p);
+	else phases(k) = jec_phase_poly(freq, p);
       }
       cout << " P4PhaseCor::readFitParamFile()/Info: Setting phases (fit-type="<<sftyp<<") for feed="<<feed<<" -> fid="<<fid<<endl;
       setFeedPhases(fid, phases);

p4phcor.h

@@ -74,6 +74,12 @@ public:
   {
     return 2.*atan(1./tan(M_PI*((p[0]*(freq+p[1]))+0.5)));
   }
+
+  /*! \brief return the phase value as a polynom function of frequency, using Jean-Eric Sawtooth fit 
+              parameters p[0,1,...] :  p[0] center freq p[1]= degree of poly, p[2].... p[n] coeffs 
+              p[2]+p[3]*(freq-p[0])+...  */
+  static double jec_phase_poly(double freq, double * p);
+  
 protected:
   // for each feed, there is a vector which contains the value of the phase as a function of frequency 
   std::vector< TVector< double > > v_phases_;

visiavg.cc

@@ -91,6 +91,8 @@ int main(int narg, const char* arg[])
   P4AVisiNumEncoder  visiencod;
   vector<sa_size_t> KVAC = visiencod.getAllAutoCor();
   vector<sa_size_t> KVCXHH = visiencod.getAllHCrossCor();
+
+  
   cout << " List of AutoCorrelation rows:"<<endl;
   for(size_t k=0; k<KVAC.size(); k++) {
     cout << "KVAC["<<k<<"]="<<KVAC[k]<<"  ->"<<visiencod.Convert2VisiName(KVAC[k])<<endl;
@@ -156,8 +158,9 @@ int main(int narg, const char* arg[])
     while (fgok) {
       //reads next visibility matrix window 
       fgok = wreader.Shift();
-      if (!fgok)  break;
       
+      if (!fgok)  break;
+     
       vismtx = wreader.getAverageVisMtx(cfdate);
       
       if (cnt==0)  {    //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations 
@@ -178,7 +181,7 @@ int main(int narg, const char* arg[])
 	//allocating 6 Cross-Corr H-H time-frequency maps 
 	cout<<"visiavg/Info: allocating H-H cross-cor Time-Frequency maps : Time->NX="<<tfmSX<<" x Freq->NY="<<tfmSY<<endl;
 	for(int k=0; k<6; k++) vtfm.push_back( TArray< complex<r_4> >(tfmSX, tfmSY) );
-	cout << "and the 2x6 for squares of rela & imag parts " << endl;
+	cout << "and the 2x6 for squares of real & imaginary parts " << endl;
 	for(int k=0; k<6; k++) vtfm_rp_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
 	for(int k=0; k<6; k++) vtfm_ip_sq.push_back( TArray< r_4 >(tfmSX, tfmSY) );
 
@@ -254,7 +257,7 @@ int main(int narg, const char* arg[])
 	
 	double tdif =  cfdate.TimeDifferenceSeconds(cfdate,dateobs)/2.;
 	timevec(TFMtmidx) = dateobs.TimeDifferenceSeconds(dateobs.ShiftSeconds (tdif ),datestart);	// centre du bin 
-	ravec(TFMtmidx) = P4Coords::RAFromTimeTU(dateobs);
+	ravec(TFMtmidx) = P4Coords::RAFromTimeTU(dateobs.ShiftSeconds (tdif ));
 	TFMtmidx++;
 	//  ... done 
 	I=0;  cntnt++;
@@ -272,7 +275,7 @@ int main(int narg, const char* arg[])
     FitsInOutFile  * fos = NULL ;
 
     if (fitsoutfile.length()>=1){
-      cout << " fitsoutfile" <<fitsoutfile<< endl;
+      cout << " fitsoutfile :" <<fitsoutfile<<":"<< endl;
       fos = new FitsInOutFile(fitsoutfile, FitsInOutFile::Fits_Create);
     }
     
@@ -350,10 +353,11 @@ int main(int narg, const char* arg[])
     //    resu.Update();
     P4FreqBand myp4fre;
     TVector <double> avg_freqs( myp4fre.getP4NbFreqChannels()/TFMfbin);
+    
     double frbase =  myp4fre.freqstart_ + myp4fre.getP4FreqResolution()/2. ;
     for (int kf=0 ; kf< myp4fre.getP4NbFreqChannels()/TFMfbin ; kf++,frbase += myp4fre.getP4FreqResolution()*TFMfbin )
       avg_freqs(kf) = frbase ;
-
+    
     potfm << PPFNameTag("FreqVec") << avg_freqs;
     if (fos != NULL) {
       ext_names.push_back("Frequences");
@@ -366,7 +370,7 @@ int main(int narg, const char* arg[])
       cout << ext_names << endl;
       delete(fos);
     }
-
+    cout << " return code "<<rc<<endl;
     cout << resu;   // Update est fait lors du print
 
   }
@@ -383,7 +387,7 @@ int main(int narg, const char* arg[])
   catch (...) {
     cerr << " visiavg.cc catched unknown (...) exception  " << endl; 
     rc = 79; 
-  } 
+  }
 
   cout << ">>>> visiavg.cc ------- END ----------- RC=" << rc << endl;
   return rc;