optimization of computation of integrated power in visi2ntac.cc, Reza 13/04/2015

visi2ntac.cc

@@ -105,6 +105,10 @@ int main(int narg, char* arg[])
     pin>>PPFNameTag("gains")>>gains;
     pin>>PPFNameTag("gn")>>gn;
     for (int p=0; p<8; p++) cout<<"gn["<<p<<"]="<<gn(p)<<" ";  
+    // on transforme gains en facteur multiplicatif
+    for(sa_size_t r=0;r<gains.NRows();r++) {
+      for(sa_size_t c=0; c<gains.NCols(); c++) gains(r,c)=1./gains(r,c);
+    }
     cout<<endl;
   }
   // Bande de frequence de 2 MHz autour du 21 cm (1420 MHz -> JF=2792) 
@@ -142,6 +146,7 @@ int main(int narg, char* arg[])
     vreader.setPrintLevel(prtlev);
     bool fgok=true;
     TMatrix< complex<r_4> > vismtx;
+    TMatrix< complex<r_4> > acsum;
     TimeStamp dateobs;
     double mttag;
     int cnt=0, cntnt=0, pcntnt=0;
@@ -149,24 +154,26 @@ int main(int narg, char* arg[])
     while (fgok) {
       fgok=vreader.ReadNext(vismtx, dateobs, mttag);
       if (fgok)  {
+	if (cnt==0)  acsum.SetSize(8, vismtx.NCols());
 	if (I==0) {
 	  xnt[0]=cnt;   xnt[1]=dateobs.SecondsPart();
 	  if (prtlev>0) 
 	    cout << "visi2ntac/Info:  dateobs="<<dateobs<<" SecondsPart()="<<dateobs.SecondsPart()<<endl;
-	  for(int k=0; k<24; k++)  xnt[2+k]=0.;
-	}
-	for(int k=0; k<8; k++)  {    
-	  TVector< complex<r_4> > vacz = vismtx.Row(KVAC[k]);
-	  Vector vac(vismtx.NCols());
-	  // correct for gain(nu)
-	  for(sa_size_t c=0; c<vac.Size(); c++) vac(c)=vacz(c).real()/gains(k,c);
-	  xnt[2+k] += vac(freqs).Sum();    // integration en bande large 
-	  xnt[2+8+k] += vac(freqs21).Sum();   // integration a +/- 1 MHz @ 1420 MHz 
-	  xnt[2+16+k] += vac(freqs21_5).Sum();   // integration a +/- 5 MHz @ 1420 MHz 
-	  xnt[2+24+k] += vac(freqsB).Sum();    // integration en bande large B
+	  for(int k=0; k<32; k++)  xnt[2+k]=0.;
+	  acsum = complex<r_4>(0.,0.);
 	}
+	for(int k=0; k<8; k++)    acsum.Row(k) += vismtx.Row(KVAC[k]);
 	I++; 
 	if (I==deltaIavg) {
+	  Vector vac(acsum.NCols());
+	  for(int k=0; k<8; k++)  {
+	  // correct for gain(nu)
+	    for(sa_size_t c=0; c<vac.Size(); c++) vac(c)=acsum(k,c).real()*gains(k,c);
+	    xnt[2+k] += vac(freqs).Sum();    // integration en bande large 
+	    xnt[2+8+k] += vac(freqs21).Sum();   // integration a +/- 1 MHz @ 1420 MHz 
+	    xnt[2+16+k] += vac(freqs21_5).Sum();   // integration a +/- 5 MHz @ 1420 MHz 
+	    xnt[2+24+k] += vac(freqsB).Sum();    // integration en bande large B
+	  }
 	  // puissance moyenne dans la bande large
 	  double fnorm=(1./(double)deltaIavg)/(double)(JFmax-JFmin+1);
 	  for(int k=0; k<8; k++)  xnt[2+k]*=fnorm;