Changeset 612 in Sophya for trunk/Poubelle

Timestamp:

Nov 22, 1999, 10:43:44 AM (26 years ago)

Author:

ansari

Message:

fin de fichier, start

Location:

trunk/Poubelle/archTOI.old

Files:

• archfileset.cc (modified) (1 diff)
• archparam.cc (modified) (1 diff)
• archvers.h (modified) (1 diff)
• gondolageom.cc (modified) (5 diffs)
• gondolageom.h (modified) (2 diffs)
• starmatcher.cc (modified) (5 diffs)
• toiboloproducer.cc (modified) (2 diffs)
• toiboloproducer.h (modified) (1 diff)
• toiiter.cc (modified) (4 diffs)

trunk/Poubelle/archTOI.old/archfileset.cc

@@ -155 +155 @@
     if (fich.nextBlock()) {
       file1stSamp[*i] = fich.blockNum()*72; // premier numsample
-      cout << "File " << *i << " 1st sample = " << fich.blockNum()*72 << endl;
+      cout << "File " << *i << " 1st sample = " << fich.blockNum()*72 
+                            << " UTC " << archParam.acq.SN2UTC(fich.blockNum()*72) << endl;
     } else {
       cout << "File " << *i << " unrecoverable, skipping" << endl;

trunk/Poubelle/archTOI.old/archparam.cc

@@ -21 +21 @@
 ArchParam::SSTParam::SSTParam() 
 :  soustPente   (true),
-   analFine     (false)
+   analFine     (true)
 {}
 

trunk/Poubelle/archTOI.old/archvers.h

@@ -2 +2 @@
 #define ARCHVERS_H
 
-#define ARCHTOI_VERS "2.0"
-#define ARCHTOI_TAG  "V_011199"
+#define ARCHTOI_VERS "3.0"
+#define ARCHTOI_TAG  "V_191199"
 
 #endif

trunk/Poubelle/archTOI.old/gondolageom.cc

@@ -4 +4 @@
 #include <math.h>
 #include "gondolageom.h"
+
+
 extern "C" {
 #include "aa_hadec.h"
@@ -27 +29 @@
             
 
-GondolaGeom::GondolaGeom() {
+GondolaGeom::GondolaGeom()
+: fit(200, 2)
+{
   azPend = 0;
   angPend = 0;
@@ -51 +55 @@
   if (nstars<0) {
     nstars = 0;
-    saz=staz=st=st2=0;
+    fit.clear();
   }
   
@@ -61 +65 @@
   if (azCor - az0 >  180) azCor -= 360;
   if (azCor - az0 < -180) azCor += 360;
-  saz  += azCor;
-  staz += azCor*deltasn;
-  st   += deltasn;
-  st2  += deltasn*deltasn;
+  
+  fit.addData(deltasn, azCor);
 }
 
@@ -70 +72 @@
 int GondolaGeom::solveStars() {
   if (nstars<2) return -1;
-  staz /= nstars;
-  st   /= nstars;
-  saz  /= nstars;
-  st2  /= nstars;
-  double a = (staz - st*saz) / (st2 - st*st);
-  double b = saz - a*st;
-  
-  nstars = -1;
-  
-  azimut = b;
+  fit.doFit();
+  azimut = fit.value(0);
   
   if (azimut > 360) azimut -= 360;
   if (azimut < 0)   azimut += 360;
 
+  nstars = -1;
+  
   return 0;
 }

trunk/Poubelle/archTOI.old/gondolageom.h

@@ -4 +4 @@
 #ifndef GONDOLAGEOM_H
 #define GONDOLAGEOM_H
+
+#include "polfitclip.h"
 
 class GondolaGeom {
@@ -59 +61 @@
   
   long nstars;
-  double saz,staz,st,st2;
   double az0;
+  
+  PolFitClip fit;
 };
 

trunk/Poubelle/archTOI.old/starmatcher.cc

@@ -11 +11 @@
 
 #define STARDUMP
+
+#define TEchan TFin
 
 #include <math.h>
@@ -271 +273 @@
       
         double da = azim-az; if (da>360) da -= 360;
-        if (da < -0.6 || da > 0.4) continue;
+       // if (da < -0.6 || da > 0.4) continue; // appropriate for TEchan
+        if (da < -0.7 || da > 0.3) continue; // appropriate for TFin
         double elv0 = elv - GondolaGeom::sstPixelHeight * lastStar.NoDiode;
         if (fabs(elv0-GondolaGeom::elevSST0) > 0.25) continue; // Might be too strong
@@ -568 +571 @@
 int StarMatcher::getPendulInfo(double sampleNum, pendulInfo& info) {
 
+  static double lastSN = -1;
+  static pendulInfo lastPendul;
+  
+  if (sampleNum == lastSN) {
+    info = lastPendul;
+    return 0;
+  }
+
   PolFitClip2 fitPendul(30,2);
 
   map<double, pendulInfo>::iterator i = pendulInfos.lower_bound(sampleNum);
   if (i == pendulInfos.begin() && (*i).second.SN >= sampleNum) return -1;
-  if (i == pendulInfos.end()   && (*i).second.SN <= sampleNum) return -1;
-  
-  if ((*i).second.SN > sampleNum) i--;
+  if (i == pendulInfos.end()) return -1;
+  map<double, pendulInfo>::iterator last = pendulInfos.end();
+  if (last == pendulInfos.begin()) return -1;
+  last--;
+  if (i == last   && (*i).second.SN <= sampleNum) return -1;
+  
+  if ((*i).second.SN > sampleNum) i--; // i just before us...
+  
+  //$CHECK$ reject if too large a gap...
+  if (sampleNum - (*i).second.SN > 1000) return -1;
+  last = i; last++;
+  if ((*last).second.SN - sampleNum > 1000) return -1;
 
   int nn=0;
   double aziprev=0, azicur=0, azi0=0;
   for (map<double, pendulInfo>::iterator ii=i; ii != pendulInfos.begin(); ii--) {
-    nn++;
     pendulInfo inf1 = (*ii).second;
+    if (fabs(inf1.SN - sampleNum) > 1000) continue;
     aziprev = azicur;
     azicur = inf1.azPendul;
+    nn++;
     if (nn==1) azi0 = azicur;
     if (nn>1 && azicur - aziprev > 180)  azicur -= 360;
@@ -593 +614 @@
   if (i != pendulInfos.end()) i++;
   for (map<double, pendulInfo>::iterator ii=i; ii != pendulInfos.end(); ii++) {
-    nn++;
     pendulInfo inf1 = (*ii).second;
+    if (fabs(inf1.SN - sampleNum) > 1000) continue;
     aziprev = azicur;
     azicur = inf1.azPendul;
+    nn++;
     if (nn>1 && azicur - aziprev > 180)  azicur -= 360;
     if (nn>1 && azicur - aziprev < -180) azicur += 360;
@@ -610 +632 @@
   if (info.azPendul <   0) info.azPendul += 360;
   info.angPendul  = fitPendul.valueY(sampleNum);
+  
+  lastSN = sampleNum;
+  lastPendul = info;
+
   return 0;
 }
-
-#if 0
-
-int StarMatcher::getPendulInfo(double sampleNum, pendulInfo& info) {
-  static double* sn    = ::vector(1, 100);
-  static double* aziP  = ::vector(1, 100);
-  static double* ampP  = ::vector(1, 100);
-  static double* sig   = ::vector(1, 100);
-  static double* aAzi  = ::vector(1, 3);
-  static double* aAmp  = ::vector(1, 3);
-  static int*   ia    = ::ivector(1,3);
-  static double** cov  = ::matrix(1,3,1,3);
-  int ndata = 0;
-  map<double, pendulInfo>::iterator i = pendulInfos.lower_bound(sampleNum);
-  if (i == pendulInfos.begin() && (*i).second.SN >= sampleNum) return -1;
-  if (i == pendulInfos.end()   && (*i).second.SN <= sampleNum) return -1;
-  
-  if ((*i).second.SN > sampleNum) i--;
-
-  int nn=0;
-  for (map<double, pendulInfo>::iterator ii=i; ii != pendulInfos.begin(); ii--) {
-    nn++;
-    ndata++;
-    pendulInfo inf1 = (*ii).second;
-    sn[ndata]   = inf1.SN;
-    ampP[ndata] = inf1.angPendul;
-    aziP[ndata] = inf1.azPendul;
-    int prev = ndata-1;
-    if (ndata>1 && aziP[ndata] - aziP[prev] > 180)  aziP[ndata] -= 360;
-    if (ndata>1 && aziP[ndata] - aziP[prev] < -180) aziP[ndata] += 360;
-    sig[ndata]  = 1;
-    if (nn>=50) break;
-  }
-  
-  nn=0;
-  if (i != pendulInfos.end()) i++;
-  for (map<double, pendulInfo>::iterator ii=i; ii != pendulInfos.end(); ii++) {
-    nn++;
-    ndata++;
-    pendulInfo inf1 = (*ii).second;
-    sn[ndata]   = inf1.SN;
-    ampP[ndata] = inf1.angPendul;
-    aziP[ndata] = inf1.azPendul;
-    int prev = ndata-1;
-    if (nn==1) prev=1;
-    if (ndata>1 && aziP[ndata] - aziP[prev] > 180)  aziP[ndata] -= 360;
-    if (ndata>1 && aziP[ndata] - aziP[prev] < -180) aziP[ndata] += 360;
-    sig[ndata]  = 1;
-    if (nn>=50) break;
-  }
-  
-  if (ndata < 3) return -1;
-  
-  ia[1] = ia[2] = ia[3] = 1;
-  double chi2;
-  try {
-    lfit(sn, aziP,  sig, ndata, aAzi, ia, 3, cov, &chi2, polfunc);
-    lfit(sn, ampP,  sig, ndata, aAmp, ia, 3, cov, &chi2, polfunc);
-  } catch(string st) {
-    return -1;
-  }
-
-  info.SN = sampleNum;
-  info.azPendul   = polval(sampleNum, aAzi, 3);
-  if (info.azPendul > 360) info.azPendul -= 360;
-  if (info.azPendul <   0) info.azPendul += 360;
-  info.angPendul  = polval(sampleNum, aAmp, 3);
-  return 0;
-}
-
-#endif
-
-#if 0
-
-int StarMatcher::getPendulInfo(double sampleNum, pendulInfo& info) {
-  map<double, pendulInfo>::iterator i = pendulInfos.lower_bound(sampleNum);
-  if (i == pendulInfos.begin() && (*i).second.SN >= sampleNum) return -1;
-  if (i == pendulInfos.end()   && (*i).second.SN <= sampleNum) return -1;
-  
-  if ((*i).second.SN > sampleNum) i--;
-  pendulInfo inf1 = (*i).second;
-  i++;
-  pendulInfo inf2 = (*i).second;
-  
-  info.SN = sampleNum;
-  if (inf2.azPendul - inf1.azPendul > 180)  inf2.azPendul -= 360;
-  if (inf2.azPendul - inf1.azPendul < -180) inf2.azPendul += 360;
-  info.azPendul  = inf1.azPendul  + (inf2.azPendul  -  inf1.azPendul) * (sampleNum - inf1.SN) / (inf2.SN - inf1.SN);
-  if (info.azPendul<0)    info.azPendul += 360;
-  if (info.azPendul>360)  info.azPendul += 360;
-  info.angPendul = inf1.angPendul + (inf2.angPendul - inf1.angPendul) * (sampleNum - inf1.SN) / (inf2.SN - inf1.SN);
-  return 0;
-}
-
-#endif
 
 

trunk/Poubelle/archTOI.old/toiboloproducer.cc

@@ -5 +5 @@
 #include "archexc.h"
 #include "requesthandler.h"
+#include "polfitclip.h"
 
-#define filt2 "boloMuV"
+#define boloMuV "boloMuV"
 
 
 TOIBoloProducer::TOIBoloProducer() {
-  possibleTOIs.insert(TOI(filt2, TOI::all, "", "microVolts"));  
+  possibleTOIs.insert(TOI(boloMuV, TOI::all, "linfilt sqfilt", "microVolts"));  
 }
 
+// No option == linfilt
+
 string TOIBoloProducer::getName() {
-  return("TOIBoloProducer 1.0");
+  return("TOIBoloProducer 1.1");
 }
+
+int TOIBoloProducer::filtHalfRange = 20;
 
 
@@ -21 +26 @@
   if (source->canGetValue(sampleNum-1, toi)) {
     value = (value + source->getValue(sampleNum-1, toi))/2.;
-    TOI toi2 = toi;
-    toi2.name = filt2;
+    
+    //TOI toi2 = toi;
+    //toi2.name = boloMuV;
+
+    TOI toi2(boloMuV, toi.index, "sqfilt", toi.unit);
     computedValue(toi2, sampleNum, value);
+  }
+  
+  // $CHECK$ possible optimization : keep fit for several successive values
+  
+  if (source->canGetValue(sampleNum-2*filtHalfRange, toi)) {
+    PolFitClip fit(2*filtHalfRange,1);
+    int s = sampleNum % 2 ? 1 : -1;
+    for (int i=0; i<2*filtHalfRange-1; i++) {
+      fit.addData(sampleNum - i, 
+         s*(source->getValue(sampleNum - i, toi) -  source->getValue(sampleNum - i -1, toi))/2);
+      s = -s;
+    }
+    fit.doFit();
+    double y = source->getValue(sampleNum - filtHalfRange, toi);
+    s = sampleNum % 2 ? 1 : -1;
+    y =y - s*fit.value(sampleNum - filtHalfRange);
+    TOI toi2(boloMuV, toi.index, "linfilt", toi.unit);
+    if (isProducing(toi2))
+      computedValue(toi2, sampleNum - filtHalfRange, y);
+    TOI toi3(boloMuV, toi.index, "", toi.unit);
+    if (isProducing(toi3))
+      computedValue(toi3, sampleNum - filtHalfRange, y);
   }
 }
 
 set<TOI> TOIBoloProducer::reqTOIFor(TOI const& toi) {
-  TOI toi2 = toi;
   set<TOI> t; 
-  if (toi.name == filt2) {
-    toi2.name = "boloRawMuV";
+  if (toi.name == boloMuV) {
+    TOI toi2("boloRawMuV", toi.index);
+    //toi2.name = "boloRawMuV";
     t.insert(toi2);
   } else {

trunk/Poubelle/archTOI.old/toiboloproducer.h

@@ -18 +18 @@
   virtual set<TOI>     reqTOIFor(TOI const&);
   virtual void         propagateLowBound(TOI const&, long sampleNum);
+  
+  static int filtHalfRange;
 };
 

trunk/Poubelle/archTOI.old/toiiter.cc

@@ -205 +205 @@
   processRequest("#UTCORIGIN 1376.5");
   processRequest("#ASIGPS ASI_GPS_archeops1999.ascii");
-  processRequest("#COMMENT Archtoi V2 -- october 1999 -- Eric Aubourg CEA/DAPNIA");
-  processRequest("#COMMENT ***WARNING***");
-  processRequest("#COMMENT ***SOME TOI'S ARE PRELIMINARY***");
-  processRequest("#COMMENT gyroSpeed is not calibrated");
-  processRequest("#COMMENT azimut/alpha/delta use galaxy crossings and ASI GPS data");
-  processRequest("#COMMENT   and assume no pendulation");
-  processRequest("#COMMENT Focal plane center elevation found at 41.5 deg");
-  processRequest("#COMMENT   with Jupiter");
-  processRequest("#COMMENT boloMuV2 is not protected against glitches");
+  processRequest("#COMMENT Archtoi V3 -- novembre 1999 -- Eric Aubourg CEA/DAPNIA");
 }
 
@@ -233 +225 @@
   }
 
-  fset.setMJDRange(mjdStart-0.01, mjdEnd+0.01);
-  fset.setSNumRange(sStart-1000, sEnd+1000);
+  // Let's add some time on each side, 30 seconds should be ok, at least
+  // for Trapani
+  
+  double delT = 30. / 86400;
+  long   delSN = long(delT / archParam.acq.perEch);
+
+  fset.setMJDRange(mjdStart-delT, mjdEnd+delT);
+  fset.setSNumRange(sStart-delSN, sEnd+delSN);
 
   fset.init();
@@ -267 +265 @@
 // Si on a epuise les fichiers de donnees, on s'arrete des qu'aucune
 // TOI n'a de valeurs apres curSample...
-  if (curSample <= 0) curSample = sStart-1;
+  if (curSample <= 0) curSample = sStart-1;  
+  if (curSample <= 0) {
+    long sn = archParam.acq.MJD2SN(mjdStart) - 1;
+    if (sn > curSample) curSample = sn;
+  }
+  
   if (curSample <= 0) curSample = fset.getSampleIndex()-1;
   
@@ -288 +291 @@
     }
     bool found=false;
-    bool valuesAhead = false;
+    //bool valuesAhead = false;
     for (vector<TOIInfo>::iterator i = request.begin(); i != request.end(); i++) {
       if (((*i).third & triggering) == 0) continue;
       if ((*i).second->canGetValue(curSample, (*i).first)) {found=true;break;}
-      if ((*i).second->lastSampleNum((*i).first)>curSample) valuesAhead=true;
+      //if ((*i).second->lastSampleNum((*i).first)>curSample) valuesAhead=true;
     }
     if (found) break;
-    if (endFound && !valuesAhead) return false;
+    if (endFound && curSample >= fset.getSampleIndex()+71) return false;
   }
   return true;