Changeset 555 in Sophya


Ignore:
Timestamp:
Nov 9, 1999, 3:04:05 PM (26 years ago)
Author:
ansari
Message:

portage cxx en cours

Location:
trunk/Poubelle/archTOI.old
Files:
1 added
26 edited

Legend:

Unmodified
Added
Removed

  • trunk/Poubelle/archTOI.old/Makefile_nodep

    r443 r555  
    22  CXX  = cxx
    33endif
    4 ifndef (CFITSIODIR)
     4ifndef CFITSIODIR
    55  CFITSIODIR = ../cfitsio
    66endif
    7 CPPFLAGS    =-O -g -I$(CFITSIODIR)
     7CPPFLAGS    =-O -g -I$(CFITSIODIR) -DANSI
    88
    99LDLIBS=-L$(CFITSIODIR)/$(HOSTTYPE) -lcfitsio
    … …  
    2121
    2222clean:
    23         rm -f archtoi *.o archdump
     23        rm -f archtoi *.o
    2424
    2525
    26 archtoi: archtoimain.o archtoi.o archeopsfile.o decompress.o toiiter.o toisvr.o archeops.o \
    27          ssthandler.o toiinterpolator.o gpsparser.o asigps.o \
     26archtoi: archtoimain.o archtoi.o archeopsfile.o archfileset.o \
     27         archeops.o decompress.o arcunit.o \
     28         tokenizer.o \
     29         toi.o toiiter.o toimanager.o toiproducer.o \
     30         toiabsorber.o toiderivproducer.o toipullproducer.o \
     31         toirepeater.o toiflagger.o toiinterpolator.o \
     32         toillboloproducer.o toilldiluproducer.o \
     33         toillgpsproducer.o toillgyroproducer.o toillreglageproducer.o \
     34         toillsstproducer.o \
     35         gpsparser.o asigps.o toiauxgpsproducer.o toiboloproducer.o \
    2836         auxinterpgps.o tsid.o archparam.o dyffttools.o fft_mayer.o \
    2937         formepulse.o pisteetoile.o sstetoile.o transfelec.o \
    30          gyrohandler.o aa_hadec.o templocator.o plgalcross.o \
    31          arcunit.o
     38         aa_hadec.o templocator.o plgalcross.o \
     39         tsidproducer.o timetoiproducer.o sststarfinder.o \
     40         starmatcher.o galcrosslocator.o gyrocalibrator.o \
     41         rotspeed.o  gondolageom.o \
     42         nrutil.o lfit.o gaussj.o covsrt.o
    3243        $(LINK.cc) $^ $(LOADLIBES) $(LDLIBS) -o $@ -lm
    33 archdump: archtoidump.o archeopsfile.o decompress.o archeops.o ssthandler.o
    34         $(LINK.cc) $^ $(LOADLIBES) $(LDLIBS) -o $@ -lm
    3544
     45

  • trunk/Poubelle/archTOI.old/archtoi.cc

    r534 r555  
    214214void ArchTOI::endLine_A() {
    215215  if (!ostr) return;
    216   *ostr << '\n';
     216  *ostr << endl;
     217 // *ostr << '\n';
    217218}
    218219

  • trunk/Poubelle/archTOI.old/archtoimain.cc

    r534 r555  
    3232    }
    3333   
    34 //   ProfilerInit(collectDetailed, bestTimeBase, 10000, 200);
     34//    ProfilerInit(collectDetailed, bestTimeBase, 10000, 200);
    3535
    3636    if (argv[1] == string("-h")) {

  • trunk/Poubelle/archTOI.old/formepulse.cc

    r534 r555  
    140140                char s2[32]="Resultatfiltrage.txt";
    141141                dlout.open(s2,ios::out | ios::trunc);
    142                 if (!dlout.is_open()){ 
    143                         cerr<<"erreur a l'ouverture du fichier de resultats: "<<s2<<endl;
    144                         exit(-1);
    145                 }
     142        //      if (!dlout.is_open()){ 
     143        //              cerr<<"erreur a l'ouverture du fichier de resultats: "<<s2<<endl;
     144        //              exit(-1);
     145        //      }
    146146
    147147                dlout<<"Impulsion courant de: "<<SSTFPulseCourant<<" Ampères"<<endl;

  • trunk/Poubelle/archTOI.old/galcrosslocator.cc

    r534 r555  
    99#define alphaZenith     "alphaZenith"
    1010#define deltaZenith     "deltaZenith"
    11 #define rotSpeed        "rotSpeed"
     11#define rotSpeed        "rotSpeed0"
    1212#define rotSpeedSample1 "rotSpeedSample1"
    1313#define rotSpeedSample2 "rotSpeedSample2"
    … …  
    2222#define alphaBolo       "alphaBolo"
    2323#define deltaBolo       "deltaBolo"
     24#define azimuthSST      "azimuthSST"
     25#define elvSST          "elvSST"
     26#define alphaSST        "alphaSST"
     27#define deltaSST        "deltaSST"
    2428
     29// Note : rotSpeed is an estimate, does not take into account balloon drift.
     30// accurate value between two galactic crossings has to be determined afterwards
     31// as a derived TOI.
    2532
    2633GalCrossLocator::GalCrossLocator() {
    27   possibleTOIs.insert(TOI(alphaZenith,      TOI::unspec, "", "hours",   "galcross0")); 
    28   possibleTOIs.insert(TOI(deltaZenith,      TOI::unspec, "", "degrees", "galcross0")); 
    29   possibleTOIs.insert(TOI(rotSpeed,         TOI::unspec, "", "deg/s",   "galcross0")); 
    30   possibleTOIs.insert(TOI(rotSpeedSample1,  TOI::unspec, "", "integer", "galcross0")); 
    31   possibleTOIs.insert(TOI(rotSpeedSample2,  TOI::unspec, "", "integer", "galcross0")); 
    32   possibleTOIs.insert(TOI(lastCrossSample,  TOI::unspec, "", "integer", "galcross0")); 
    33   possibleTOIs.insert(TOI(lastCrossSample,  TOI::unspec, "", "integer", "galcross0")); 
    34   possibleTOIs.insert(TOI(nextCrossSample,  TOI::unspec, "", "integer", "galcross0")); 
    35   possibleTOIs.insert(TOI(azimuthFPC,       TOI::unspec, "", "degrees", "galcross0")); 
    36   possibleTOIs.insert(TOI(elvFPC,           TOI::unspec, "", "degrees", "galcross0")); 
    37   possibleTOIs.insert(TOI(alphaFPC,         TOI::unspec, "", "hours",   "galcross0")); 
    38   possibleTOIs.insert(TOI(deltaFPC,         TOI::unspec, "", "degrees", "galcross0")); 
    39   possibleTOIs.insert(TOI(azimuthBolo,      TOI::all,    "", "degrees", "galcross0")); 
    40   possibleTOIs.insert(TOI(elvBolo,          TOI::all,    "", "degrees", "galcross0")); 
    41   possibleTOIs.insert(TOI(alphaBolo,        TOI::all,    "", "hours",   "galcross0")); 
    42   possibleTOIs.insert(TOI(deltaBolo,        TOI::all,    "", "degrees", "galcross0")); 
     34  possibleTOIs.insert(TOI(alphaZenith,      TOI::unspec, "interp flag", "hours",   "")); 
     35  possibleTOIs.insert(TOI(deltaZenith,      TOI::unspec, "interp flag", "degrees", "")); 
     36  possibleTOIs.insert(TOI(rotSpeed,         TOI::unspec, "interp flag", "deg/s",   "galcross0")); 
     37  possibleTOIs.insert(TOI(rotSpeedSample1,  TOI::unspec, "interp flag", "integer", "galcross0")); 
     38  possibleTOIs.insert(TOI(rotSpeedSample2,  TOI::unspec, "interp flag", "integer", "galcross0")); 
     39  possibleTOIs.insert(TOI(lastCrossSample,  TOI::unspec, "interp flag", "integer", "galcross0")); 
     40  possibleTOIs.insert(TOI(nextCrossSample,  TOI::unspec, "interp flag", "integer", "galcross0")); 
     41  possibleTOIs.insert(TOI(azimuthFPC,       TOI::unspec, "interp flag", "degrees", "galcross0")); 
     42  possibleTOIs.insert(TOI(elvFPC,           TOI::unspec, "interp flag", "degrees", "galcross0")); 
     43  possibleTOIs.insert(TOI(alphaFPC,         TOI::unspec, "interp flag", "hours",   "galcross0")); 
     44  possibleTOIs.insert(TOI(deltaFPC,         TOI::unspec, "interp flag", "degrees", "galcross0")); 
     45  possibleTOIs.insert(TOI(azimuthBolo,      TOI::all,    "interp flag", "degrees", "galcross0")); 
     46  possibleTOIs.insert(TOI(elvBolo,          TOI::all,    "interp flag", "degrees", "galcross0")); 
     47  possibleTOIs.insert(TOI(alphaBolo,        TOI::all,    "interp flag", "hours",   "galcross0")); 
     48  possibleTOIs.insert(TOI(deltaBolo,        TOI::all,    "interp flag", "degrees", "galcross0")); 
     49  possibleTOIs.insert(TOI(azimuthSST,       TOI::unspec, "interp flag", "degrees", "galcross0")); 
     50  possibleTOIs.insert(TOI(elvSST,           TOI::unspec, "interp flag", "degrees", "galcross0")); 
     51  possibleTOIs.insert(TOI(alphaSST,         TOI::unspec, "interp flag", "hours",   "galcross0")); 
     52  possibleTOIs.insert(TOI(deltaSST,         TOI::unspec, "interp flag", "degrees", "galcross0")); 
    4353}
    4454
    … …  
    5565double GalCrossLocator::getValue(long sampleNum, TOI const& toi) {
    5666  if (!canGetValue(sampleNum, toi)) return -1;
     67 
     68  if (toi.options.find("flag") != toi.options.end()) {
     69     int SN1, SN2;
     70     tempLocator.getCrossSamples(sampleNum, SN1, SN2);
     71     if (sampleNum == SN1 || sampleNum == SN2) return 1;
     72     return 0;
     73  }
    5774
    58   map<TOI, TOIProducer*> m = neededTOIs[toi];
     75  map<TOI, TOIProducer*>  & m = neededTOIs[toi];
    5976  double lat, lon, ts;
    6077 
    6178  for (map<TOI, TOIProducer*>::iterator i = m.begin(); i != m.end(); i++) {
    62     TOI inToi = (*i).first;
     79    TOI const& inToi = (*i).first;
    6380    TOIProducer* prod =  (*i).second;
    6481    if (inToi.name == "latitude")  lat = prod->getValue(sampleNum, inToi);
    … …  
    83100  if (toi.name == alphaBolo)       return tempLocator.getAlphaBolo(sampleNum, toi.index);
    84101  if (toi.name == deltaBolo)       return tempLocator.getDeltaBolo(sampleNum, toi.index);
     102
     103  if (toi.name == azimuthSST)      return tempLocator.getAzimutSST(sampleNum);
     104  if (toi.name == elvSST)          return tempLocator.getElvSST(sampleNum);
     105  if (toi.name == alphaSST)        return tempLocator.getAlphaSST(sampleNum);
     106  if (toi.name == deltaSST)        return tempLocator.getDeltaSST(sampleNum);
    85107
    86108  if (toi.name == lastCrossSample || toi.name == rotSpeedSample1) {

  • trunk/Poubelle/archTOI.old/gyrocalibrator.cc

    r534 r555  
    44// assumption : same calibration for all gyros...
    55
     6// Warning, current implementation can only output ONE calibration.
     7// workaround : clone the object for different options.
     8
    69
    710#include "gyrocalibrator.h"
     11#include "archexc.h"
     12#include "archparam.h"
    813
    914#define gyroCal    "gyroCal"
     15#define gyroOffset "gyroOffset"
    1016#define gyroSpeed  "gyroSpeed"
    1117
    1218GyroCalibrator::GyroCalibrator() {
    13   possibleTOIs.insert(TOI(gyroCal,      TOI::unspec, "linearCal", "deg/s/V")); 
     19  possibleTOIs.insert(TOI(gyroCal,      TOI::all,    "linearCal", "deg/s/V")); 
     20  possibleTOIs.insert(TOI(gyroOffset,   TOI::all,    "linearCal", "deg/s/V")); 
    1421  possibleTOIs.insert(TOI(gyroSpeed,    TOI::all   , "linearCal", "deg/s")); 
    1522 
    16   //needAfter = 3600;  // about one turn...
     23  needAfter = 3600;  // about one turn...
    1724  startSample = -1;
     25  lastRotSpeed = -1;
    1826  lastFence1 = lastFence2 = -1;
     27  for (int i=0; i<3; i++) {
     28    lastCalib[i] = -99999;
     29    lastOffset[i] = 0;
     30  }
     31  gyroProducer = NULL;
    1932}
    2033
    … …  
    3548}
    3649
    37 void GyroCalibrator::dataFeed(TOIProducer* source, TOI const& toi, long sampleNum, double value) {
    38   if (toi.name != "gyroV" || toi.index != 2) return;
     50void GyroCalibrator::dataFeed(TOIProducer* , TOI const& toi, long sampleNum, double value) {
     51 // if (toi.name != "gyroV" || toi.index != 2) return;
     52  if (toi.name != "gyroV" ) return;
    3953 
    40   if (gyro2.empty()) {
     54  if (gyro[2].empty()) {
    4155    startSample = sampleNum;
    4256  }
    4357 
    44   gyro2.push_back(value);
     58  gyro[toi.index].push_back(value);
     59}
     60
     61void GyroCalibrator::recomputeCalib() {
     62  // Integral of gyro signal between the two fences.
     63 
     64  // can we do it ?
     65  for (int i=0; i<3; i++) {
     66    lastCalib[i] = -99999;
     67  }
     68  if (startSample > lastFence1) return;
     69  if ((long)gyro[2].size()+startSample < lastFence2) return;
     70       
     71  double sum[3];
     72  for (int i=0; i<3; i++) {
     73    sum[i] = 0;
     74  }
     75  vector<double>::iterator i0     = gyro[0].begin() + (lastFence1-startSample);
     76  vector<double>::iterator i1     = gyro[1].begin() + (lastFence1-startSample);
     77  vector<double>::iterator i2     = gyro[2].begin() + (lastFence1-startSample);
     78  vector<double>::iterator stop2  = gyro[2].begin() + (lastFence2-startSample);
     79  for (; i2 != stop2; i0++,i1++,i2++) {
     80    sum[0] += *i0;
     81    sum[1] += *i1;
     82    sum[2] += *i2;
     83  }
     84   
     85  // average of signal
     86  sum[0] /= (lastFence2-lastFence1);
     87  sum[1] /= (lastFence2-lastFence1);
     88  sum[2] /= (lastFence2-lastFence1);
     89 
     90  lastOffset[0] = sum[0];
     91  lastOffset[1] = sum[1];
     92  lastOffset[2] = (lastOffset[0] + lastOffset[1])/2.;  // We don't have a better estimate...
     93 
     94  lastCalib[2]  = lastRotSpeed / (sum[2] - lastOffset[2]);
     95  lastCalib[0] = lastCalib[1] = lastCalib[2];          // We don't have a better estimate...
     96}
     97
     98bool GyroCalibrator::fetchFences(long sampleNum) {
     99  map<TOI, TOIProducer*> & m = (*(neededTOIs.begin())).second;
     100 
     101  long oldf2 = lastFence2;
     102  for (map<TOI, TOIProducer*>::iterator i = m.begin(); i != m.end(); i++) {
     103    TOI const& inToi = (*i).first;
     104    TOIProducer* prod =  (*i).second;
     105    if (inToi.name == "rotSpeed")        lastRotSpeed = prod->getValue(sampleNum, inToi);
     106    if (inToi.name == "rotSpeedSample1") lastFence1   = prod->getValue(sampleNum, inToi);
     107    if (inToi.name == "rotSpeedSample2") lastFence2   = prod->getValue(sampleNum, inToi);
     108    if (inToi.name == "gyroV")           gyroProducer = prod;
     109  }
     110 
     111  return (oldf2 != lastFence2);
    45112}
    46113
    47114
     115bool GyroCalibrator::canGetValue(long sampleNum, TOI const& ) {
     116  // We can get value if sampleNum is between the two fences, or
     117  // if a new fence later than sampleNum is now available.
     118  // In any case, we must have gyro data up to the highest fence.
     119 
     120  if (startSample > lastFence1) return false; // will never work...
     121  if (sampleNum >= lastFence1 && sampleNum < lastFence2 && (gyro[2].size()+startSample >= lastFence2)) return true;
     122  if ((long)gyro[2].size()+startSample < lastFence2) return false; // We have to wait for more gyro data
     123  if (fetchFences(sampleNum)) {
     124    recomputeCalib(); // do it now to keep a consistent state
     125  }
     126  if (sampleNum >= lastFence1 && sampleNum < lastFence2) {
     127    return true;
     128  } else {
     129    return false;
     130  } 
     131}
    48132
     133bool GyroCalibrator::canGetValueLater(long sampleNum, TOI const& ) {
     134  if (sampleNum >= lastFence1 && sampleNum < lastFence2 && ((long)gyro[2].size()+startSample >= lastFence2)) return false;
     135      // because can get now
     136  if (sampleNum >= lastFence2) { // check if new fence...
     137    if (fetchFences(sampleNum)) {
     138      recomputeCalib(); // do it now to keep a consistent state
     139    }
     140  }
     141  if (lastFence1<0 || lastFence2<0) return true;
     142  if (sampleNum > lastFence2) return true;
     143  return (startSample <= lastFence1 && (long)gyro[2].size()+startSample < lastFence2 && sampleNum > lastFence1);
     144}
     145
     146double GyroCalibrator::getValue(long sampleNum, TOI const& toi) {
     147  if (startSample > lastFence1) return -1; // will never work...
     148  if ((long)gyro[2].size()+startSample < lastFence2) return -1; // We have to wait for more gyro data
     149  if (!(sampleNum >= lastFence1 && sampleNum < lastFence2)) {
     150    if (fetchFences(sampleNum)) {
     151      if ((long)gyro[2].size()+startSample < lastFence2) return -1; // We have to wait for more gyro data
     152      recomputeCalib(); // do it now to keep a consistent state
     153    }
     154  }
     155 
     156  if (!(sampleNum >= lastFence1 && sampleNum < lastFence2))
     157    return -1;
     158   
     159  if (lastCalib[0] < 0) recomputeCalib();
     160   
     161  if (toi.name == gyroCal)    return lastCalib[toi.index];
     162  if (toi.name == gyroOffset) return lastOffset[toi.index];
     163  if (toi.name == gyroSpeed) {
     164    TOI toi2 = TOI("gyroV", toi.index); // $CHECK$ maybe should get from reqtoi ?
     165    double gv = gyroProducer->getValue(sampleNum, toi2);
     166    return (gv-lastOffset[toi.index])*lastCalib[toi.index];
     167  }
     168  throw ArchExc("Cannot produce "+toi.fullName());
     169}
     170
     171
     172void GyroCalibrator::propagateLowBound(TOI const&, long /*sampleNum*/) {
     173   if (startSample < lastFence1) {
     174     if (gyro[0].size() > lastFence1 - startSample) {
     175       for (int i=0; i<3; i++) {
     176         vector<double>::iterator x = gyro[i].begin() + lastFence1 - startSample;
     177         gyro[i].erase(gyro[i].begin(), x);
     178       }
     179       startSample = lastFence1;
     180     }
     181   }
     182}
     183

  • trunk/Poubelle/archTOI.old/gyrocalibrator.h

    r534 r555  
    55#define GYROCALIBRATOR_H
    66
    7 #include "toiderivproducer.h"
     7#include "toipullproducer.h"
    88
    99#include <vector>
    1010
    11 class GyroCalibrator : public TOIDerivProducer {
     11class GyroCalibrator : public TOIPullProducer {
    1212public:
    1313  GyroCalibrator();
    … …  
    1616
    1717  virtual bool   canGetValue(long sampleNum, TOI const& toi); // for this samplenum
     18  virtual bool   canGetValueLater(long sampleNum, TOI const& toi);  // Might, later but not now
     19 
    1820  virtual double getValue(long sampleNum, TOI const& toi);
    19 
     21 
    2022protected:
    2123  virtual set<TOI>     reqTOIFor(TOI const&);
    2224  virtual void         propagateLowBound(TOI const&, long sampleNum);
    2325 
    24 private:
     26  virtual void recomputeCalib();
     27  virtual bool fetchFences(long sampleNum);
     28 
     29 
    2530  long startSample;
    26   vector<double> gyro2;
     31  vector<double> gyro[3];
    2732  long lastFence1, lastFence2;
     33  double lastRotSpeed; // on gyro 2
     34  double lastCalib[3];
     35  double lastOffset[3];
     36  TOIProducer* gyroProducer;
    2837};
    2938

  • trunk/Poubelle/archTOI.old/mkd

    r430 r555  
    22
    33cp Makefile_nodep Makefile
     4if( ! $?CFITSIODIR ) then
     5  setenv CFITSIODIR ../cfitsio
     6endif
    47foreach f (*.c *.cc)
    5   gcc -MM -I../cfitsio $f >> Makefile
     8  gcc -MM -I$CFITSIODIR $f >> Makefile
    69end

  • trunk/Poubelle/archTOI.old/pisteetoile.cc

    r534 r555  
    1 #include "math.h"
     1#include <math.h>
    22#include <new.h>
    3 #include "iostream.h"
     3#include <iostream.h>
    44#include "archparam.h"
    55#include "pisteetoile.h"

  • trunk/Poubelle/archTOI.old/sststarfinder.cc

    r534 r555  
    1515#define sstStarSN    "sstStarSN"
    1616#define sstStarUTC   "sstStarUTC"
     17
     18bool  SSTStarFinder::has2bars=false;
     19int   SSTStarFinder::elecOffset=0;
     20
    1721
    1822
    … …  
    9296
    9397void SSTStarFinder::dataFeed(long sampleNum, int* diodeSignal) {
    94   if (producedTOIs.empty()) return;
     98  if (producedTOIs.empty() && processors.empty()) return;
    9599  if (sampleNum == lastFedSample+1) {
    96100        for(int i=0;i<NbPhotDiodBarette; i++) {

  • trunk/Poubelle/archTOI.old/starmatcher.cc

    r534 r555  
     1// starmatcher.cc
     2// Eric Aubourg         CEA/DAPNIA/SPP   novembre 1999
     3
    14#include "starmatcher.h"
    25#include "sststarfinder.h"
    36#include "toimanager.h"
    47#include "archexc.h"
     8#include "archparam.h"
     9#include "gondolageom.h"
    510
    611extern "C" {
    712#include "aa_hadec.h"
     13#define NRANSI
     14#include "nrutil.h"
     15
     16void lfit(float x[], float y[], float sig[], int ndat, float a[], int ia[],
     17        int ma, float **covar, float *chisq, void (*funcs)(float, float [], int));
     18
     19void polfunc(float x, float afunc[], int ma);
     20void sinfunc(float x, float afunc[], int ma);
     21}
     22
     23void polfunc(float x, float afunc[], int ma) {
     24  afunc[1] = 1;
     25  for (int i=2; i<=ma; i++)
     26    afunc[i] = afunc[i-1]*x;
     27}
     28
     29void sinfunc(float x, float afunc[], int /*ma*/) {
     30  afunc[1] = cos(x);
     31  afunc[2] = sin(x);
     32  afunc[3] = 1;
     33}
     34
     35
     36float polval(float x, float a[], int ma);
     37
     38float polval(float x, float a[], int ma) {
     39  float r = a[ma];
     40  for (int i=ma-1; i>0; i--) {
     41    r = r*x+a[i];
     42  }
     43  return r;
    844}
    945
    1046#include <stdio.h>
    1147
     48#ifdef __DECCXX
     49#define SWAP
     50#endif
     51#if defined(Linux)  || defined(linux)
     52#define SWAP
     53#endif
     54
     55typedef unsigned int4 uint_4;
     56typedef unsigned short uint_2;
     57
     58static inline void bswap4(void* p)
     59{
     60  uint_4 tmp = *(uint_4*)p;
     61  *(uint_4*)p = ((tmp >> 24) & 0x000000FF) |
     62                ((tmp >> 8)  & 0x0000FF00) |
     63                ((tmp & 0x0000FF00) << 8)  |
     64                ((tmp & 0x000000FF) << 24);
     65}
     66
     67static inline void bswap2(void* p)
     68{
     69  uint_2 tmp = *(uint_2*)p;
     70  *(uint_2*)p = ((tmp >> 8) & 0x00FF) |
     71                ((tmp & 0x00FF) << 8);
     72}
     73
     74
     75#define azimuthPendul   "azimuthPendul"
     76#define anglePendul     "anglePendul"
     77#define azimuthAxis     "azimuthAxis"
     78#define elvAxis         "deltaZenith"
     79#define alphaAxis       "alphaZenith"
     80#define deltaAxis       "deltaZenith"
     81#define azimuthFPC      "azimuthFPC"
     82#define elvFPC          "elvFPC"
     83#define alphaFPC        "alphaFPC"
     84#define deltaFPC        "deltaFPC"
     85#define azimuthBolo     "azimuthBolo"
     86#define elvBolo         "elvBolo"
     87#define alphaBolo       "alphaBolo"
     88#define deltaBolo       "deltaBolo"
     89#define azimuthSST      "azimuthSST"
     90#define elvSST          "elvSST"
     91#define alphaSST        "alphaSST"
     92#define deltaSST        "deltaSST"
     93
     94
    1295StarMatcher::StarMatcher() {
    13   possibleTOIs.insert(TOI("dummatcher",      TOI::unspec, "", "blabla")); 
    14  
     96  possibleTOIs.insert(TOI(azimuthSST,       TOI::unspec, "interp", "degrees","sstmatch")); 
     97  possibleTOIs.insert(TOI(elvSST,           TOI::unspec, "interp", "degrees","sstmatch")); 
     98  possibleTOIs.insert(TOI(alphaSST,         TOI::unspec, "interp", "hours","sstmatch")); 
     99  possibleTOIs.insert(TOI(deltaSST,         TOI::unspec, "interp", "degrees","sstmatch")); 
     100  possibleTOIs.insert(TOI(azimuthAxis,      TOI::unspec, "interp", "degrees","sstmatch")); 
     101  possibleTOIs.insert(TOI(elvAxis,          TOI::unspec, "interp", "degrees","sstmatch")); 
     102  possibleTOIs.insert(TOI(alphaAxis,        TOI::unspec, "interp", "hours","sstmatch")); 
     103  possibleTOIs.insert(TOI(deltaAxis,        TOI::unspec, "interp", "degrees","sstmatch")); 
     104  possibleTOIs.insert(TOI(azimuthPendul,    TOI::unspec, "interp", "degrees","sstmatch")); 
     105  possibleTOIs.insert(TOI(anglePendul,      TOI::unspec, "interp", "degrees","sstmatch")); 
     106  possibleTOIs.insert(TOI(azimuthFPC,       TOI::unspec, "interp", "degrees", "sstmatch")); 
     107  possibleTOIs.insert(TOI(elvFPC,           TOI::unspec, "interp", "degrees", "sstmatch")); 
     108  possibleTOIs.insert(TOI(alphaFPC,         TOI::unspec, "interp", "hours",   "sstmatch")); 
     109  possibleTOIs.insert(TOI(deltaFPC,         TOI::unspec, "interp", "degrees", "sstmatch")); 
     110  possibleTOIs.insert(TOI(azimuthBolo,      TOI::all,    "interp", "degrees", "sstmatch")); 
     111  possibleTOIs.insert(TOI(elvBolo,          TOI::all,    "interp", "degrees", "sstmatch")); 
     112  possibleTOIs.insert(TOI(alphaBolo,        TOI::all,    "interp", "hours",   "sstmatch")); 
     113  possibleTOIs.insert(TOI(deltaBolo,        TOI::all,    "interp", "degrees", "sstmatch")); 
     114
    15115  FILE* f;
    16116 
    … …  
    18118  if (!f) throw ArchExc("Error opening gsc7.dat");
    19119 
    20   fread(&nstars, sizeof(long), 1 , f);
     120  int4 n4;
     121  fread(&n4, sizeof(int4), 1 , f);
     122
     123#ifdef SWAP
     124  bswap4(&n4);
     125#endif
     126  nstars = n4;
     127
    21128  stars = new gscStar[nstars];
    22   fread(stars, sizeof(gscStar), nstars, f);
     129  char* compdata = new char[10*nstars];
     130  fread(compdata, 10, nstars, f);
    23131  fclose(f);
     132 
     133  for (int i=0; i<nstars; i++) {
     134#ifdef SWAP
     135     ((char*)&(stars[i].ra))[0]  = compdata[10*i+3];
     136     ((char*)&(stars[i].ra))[1]  = compdata[10*i+2];
     137     ((char*)&(stars[i].ra))[2]  = compdata[10*i+1];
     138     ((char*)&(stars[i].ra))[3]  = compdata[10*i+0];
     139     ((char*)&(stars[i].dec))[0] = compdata[10*i+7];
     140     ((char*)&(stars[i].dec))[1] = compdata[10*i+6];
     141     ((char*)&(stars[i].dec))[2] = compdata[10*i+5];
     142     ((char*)&(stars[i].dec))[3] = compdata[10*i+4];
     143     ((char*)&(stars[i].mag))[0] = compdata[10*i+9];
     144     ((char*)&(stars[i].mag))[1] = compdata[10*i+8];
     145#else
     146     ((char*)&(stars[i].ra))[0]  = compdata[10*i+0];
     147     ((char*)&(stars[i].ra))[1]  = compdata[10*i+1];
     148     ((char*)&(stars[i].ra))[2]  = compdata[10*i+2];
     149     ((char*)&(stars[i].ra))[3]  = compdata[10*i+3];
     150     ((char*)&(stars[i].dec))[0] = compdata[10*i+4];
     151     ((char*)&(stars[i].dec))[1] = compdata[10*i+5];
     152     ((char*)&(stars[i].dec))[2] = compdata[10*i+6];
     153     ((char*)&(stars[i].dec))[3] = compdata[10*i+7];
     154     ((char*)&(stars[i].mag))[0] = compdata[10*i+8];
     155     ((char*)&(stars[i].mag))[1] = compdata[10*i+9];
     156#endif
     157  }
     158 
     159  delete[] compdata;
    24160
    25161  TOIProducer* prod = TOIManager::findTOIProducer(TOI("sstStarCount"));
    … …  
    35171  }
    36172 
     173  lastSeq = 0;
     174 
    37175  sprod->registerProcessor(this);
    38176
    … …  
    43181}
    44182
     183#ifdef STARDUMP
    45184static ofstream starstream("stars.dat");
     185static ofstream cstarstream("cstars.dat");
     186static ofstream pendstream("pendul.dat");
     187#endif
    46188
    47189void StarMatcher::dataFeed(SSTEtoile const& x) {
    48   lastStars[(long)x.TEchan] = x;
     190  lastStars.push_back(x);
     191}
     192
     193static long lastCleanSave=0;
     194
     195void nrerror(char * error_text) {
     196  throw(string(error_text));
     197}
     198
     199
     200void StarMatcher::processStars() {
     201
     202  if (lastStars.empty()) return;
     203
     204  map<TOI, TOIProducer*> & m = (*neededTOIs.begin()).second;
     205  while (!lastStars.empty()) {
     206    SSTEtoile lastStar = lastStars.front();
     207    lastStars.pop_front();
     208
     209    double lat, lon, ts, alpha, delta, az, rspeed;
     210   
     211    long snstar = (long) lastStar.TEchan;
     212 
     213    for (map<TOI, TOIProducer*>::iterator i = m.begin(); i != m.end(); i++) {
     214      TOI const& inToi = (*i).first;
     215      TOIProducer* prod =  (*i).second;
     216      if (inToi.name == "latitude")    lat = prod->getValue(snstar, inToi);
     217      if (inToi.name == "longitude")   lon = prod->getValue(snstar, inToi);
     218      if (inToi.name == "tsid")         ts = prod->getValue(snstar, inToi);
     219      if (inToi.name == "alphaSST")  alpha = prod->getValue(snstar, inToi);
     220      if (inToi.name == "deltaSST")  delta = prod->getValue(snstar, inToi);
     221      if (inToi.name == "azimuthSST")   az = prod->getValue(snstar, inToi);
     222      if (inToi.name == "rotSpeed") rspeed = prod->getValue(snstar, inToi);
     223    }
     224   
     225    // correct azimuth using fractional value of TEchan
     226   
     227    az -= rspeed * archParam.acq.perEch * (lastStar.TEchan-snstar);
     228   
     229    // find all stars +- 2 deg boresight
     230    double dist = 2;
     231    double dmin = delta - dist; if (dmin<-90) dmin=-90;
     232    double dmax = delta + dist; if (dmax> 90) dmax= 90;
     233    double amin = alpha - dist / cos(delta * 3.1415926/180) / 15.;
     234    if (amin<0) amin += 24;
     235    double amax = alpha + dist / cos(delta * 3.1415926/180) / 15.;
     236    if (amax>24) amax -= 24;
     237 
     238    int a,b,c;
     239    a=0; c=nstars-1;
     240    while (a+1<c) {
     241      b = (a+c)/2;
     242      if (stars[b].dec < dmin) a=b; else c=b;
     243    }
     244    int imin = a;
     245    a=0; c=nstars;
     246    while (a+1<c) {
     247      b = (a+c)/2;
     248      if (stars[b].dec < dmax) a=b; else c=b;
     249    }
     250    int imax = c;
     251
     252    for (int i=imin; i<=imax; i++) {
     253      if (stars[i].ra >= amin && stars[i].ra <= amax) {
     254        double ha = (ts/3600. - stars[i].ra) * 15. * 3.1415926/180.;
     255        double elv, azim;
     256        hadec_aa(lat * 3.1415926/180., ha, stars[i].dec * 3.1415926/180.,
     257                 &elv, &azim);
     258        elv  *= 180/3.1415926;
     259        azim *= 180/3.1415926;
     260        if (azim<0) azim += 360;
     261     
     262        double da = azim-az; if (da>360) da -= 360;
     263        if (da < -0.6 || da > 0.4) continue;
     264        double elv0 = elv - 1.41/45. * lastStar.NoDiode;
     265        if (fabs(elv0-GondolaGeom::elevSST0) > 0.25) continue; // Might be too strong
     266       
     267#ifdef STARDUMP     
     268        starstream << setprecision(10)  << lastStar.TEchan << " " <<
     269          lastStar.NoDiode << " " <<
     270          alpha << " " << delta << " " <<
     271          az << " " <<
     272          stars[i].ra << " " << stars[i].dec << " " <<
     273          elv << " " << azim << " " <<
     274          lastStar.InpCurrent << " " << stars[i].mag << "\n";
     275#endif
     276         
     277        matchStar s;
     278        lastSeq++;
     279        s.SN       = lastStar.TEchan;
     280        s.raGSC    = stars[i].ra;
     281        s.decGSC   = stars[i].dec;
     282        s.azGSC    = azim;
     283        s.elvGSC   = elv;
     284        s.nDiode   = lastStar.NoDiode;
     285        s.ok       = true;
     286        s.seq      = lastSeq;
     287        s.lon      = lon;
     288        s.lat      = lat;
     289        s.ts       = ts;
     290       
     291        matchStars.push_back(s);
     292      }
     293    }
     294  }
     295 
     296  // new set of matched stars... Clean, and get parameters...
     297  // We don't want more than 20 seconds of data
     298 
     299  if (matchStars.empty()) return;
     300 
     301 
     302  double snEnd = matchStars.back().SN;
     303  deque<matchStar>::iterator it;
     304  for (it = matchStars.begin();  it!=matchStars.end(); it++) {
     305     if ((snEnd - (*it).SN)*archParam.acq.perEch < 20)
     306       break;
     307  }
     308  if (it != matchStars.begin()) {
     309    matchStars.erase(matchStars.begin(), it);
     310  }
     311 
     312  // we want to clean on the last 5 seconds of data.
     313 
     314  int nskip=0;
     315  for (it = matchStars.begin();  it!=matchStars.end(); it++) {
     316     if ((snEnd - (*it).SN)*archParam.acq.perEch < 7)
     317       break;
     318     nskip++;
     319  }
     320 
     321  if (matchStars.size()-nskip < 30) return; // pas assez d'etoiles
     322 
     323  // we remove "bursts" of stars, ie more than 4 stars in the same samplenum
     324 
     325  long lastSN = 0;
     326  deque<matchStar>::iterator lastIt = it;
     327  long burstLen = 0;
     328  for (deque<matchStar>::iterator it1 = it ;  it1!=matchStars.end(); it1++) {
     329    matchStar s = (*it1);
     330    if ((long) s.SN == lastSN) {
     331      burstLen++;
     332      continue;
     333    }
     334    if (burstLen >= 4) {
     335      for (deque<matchStar>::iterator it2=lastIt; it2 != it1; it2++) {
     336        (*it2).ok=false;
     337      }
     338    }
     339    burstLen = 1;
     340    lastIt = it1;
     341    lastSN = s.SN;
     342  }
     343  // we fit the data to a polynomial, with clipping...
     344 
     345  float* sn   =  ::vector(1, matchStars.size());
     346  float* elv0 =  ::vector(1, matchStars.size());
     347  float* azi  =  ::vector(1, matchStars.size());
     348  float* sig  =  ::vector(1, matchStars.size());
     349  float* ae   =  ::vector(1, 3);
     350  float* aa   =  ::vector(1, 3);
     351  int*  ia    = ivector(1, 3);
     352  float** cov =  matrix(1, 3, 1, 3);
     353  int ndata;
     354 
     355  long sn0 = matchStars.front().SN;
     356  long snmin;
     357  long snmax;
     358  for (int i=0; i<4; i++) {
     359    ndata = 0;
     360    snmin = 99999999;
     361    snmax = -99999999;
     362    for (deque<matchStar>::iterator it1 = it ;  it1!=matchStars.end(); it1++) {
     363      matchStar s = (*it1);
     364      if (!s.ok) continue;
     365      double delv, daz;
     366      if (i) {
     367        delv = polval(s.SN-sn0, ae, 3)-(s.elvGSC - s.nDiode*1.41/45.);
     368        daz  = polval(s.SN-sn0, aa, 3)- s.azGSC;
     369        if (daz>=180) daz -= 360;
     370        if (daz<-180) daz += 360;
     371      }
     372      double dcutelv=0.2;
     373      double dcutaz =0.4;
     374      if (i>=2) {
     375         dcutelv = 0.05;
     376         dcutaz  = 0.1;
     377      }
     378      if (i>=3) {
     379         dcutelv = 0.02;
     380         dcutaz  = 0.03;
     381      }
     382      if (i == 0 || ((fabs(delv)<dcutelv && fabs(daz)<dcutaz) && i<3)) {
     383          ndata++;
     384          sn  [ndata] = s.SN-sn0;
     385          elv0[ndata] = s.elvGSC - s.nDiode*1.41/45.;
     386          azi [ndata] = s.azGSC; // $CHECK$  ou ajuster difference avec az galcross ?
     387          if (ndata>1 && azi[ndata] - azi[ndata-1] > 180) azi[ndata] -= 360;
     388          if (ndata>1 && azi[ndata] - azi[ndata-1] < -180) azi[ndata] += 360;
     389          sig [ndata] = 0.01;
     390          if (s.SN-sn0 > snmax) snmax = s.SN-sn0;
     391          if (s.SN-sn0 < snmin) snmin = s.SN-sn0;
     392      }
     393      if ((fabs(delv)>=dcutelv || fabs(daz)>=dcutaz) && i==3) {
     394          (*it1).ok = false;
     395      }
     396    }
     397    if (i==3) break;
     398    if (ndata<5) return; // on ne peut rien faire
     399    ia[1]   =  ia[2] = ia[3] = 1;
     400    float chi2;
     401    try{
     402      lfit(sn, elv0, sig, ndata, ae, ia, 3, cov, &chi2, polfunc);
     403      lfit(sn, azi,  sig, ndata, aa, ia, 3, cov, &chi2, polfunc);
     404    } catch(string s) {
     405      return;
     406    }
     407  }
     408 
     409  for (deque<matchStar>::iterator it1 = it ;  it1!=matchStars.end(); it1++) {
     410    if ((*it1).ok && (*it1).seq > lastCleanSave) {
     411      lastCleanSave = (*it1).seq;
     412#ifdef STARDUMP
     413      cstarstream <<  (*it1).seq << "\n";
     414#endif
     415      posInfo info;
     416      info.SN      =   (*it1).SN;
     417      info.azStar  =   (*it1).azGSC;
     418      info.elvStar =   (*it1).elvGSC;
     419      info.diodStar=   (*it1).nDiode;
     420      info.lon     =   (*it1).lon;
     421      info.lat     =   (*it1).lat;
     422      info.ts      =   (*it1).ts;
     423      posInfos[info.SN] = info;
     424    }
     425  }
     426 
     427  // On a des etoiles nettoyees, on va trouver amplitude et phase du
     428  // signal en elevation, ce qui va nous donner les deux angles d'Euler
     429  // de la pendulation (au premier ordre en theta)
     430 
     431  // Il faut avoir une periode entiere ou pas loin, sinon on ne peut
     432  // rien dire simplement....
     433 
     434  it = matchStars.end(); it--;
     435  if ((((*it).SN) - (*matchStars.begin()).SN)*archParam.acq.perEch < 17) return;
     436 
     437  long snmid = (((*it).SN) + (*matchStars.begin()).SN)/2;
     438   
     439  ndata=0; 
     440 
     441  for (deque<matchStar>::iterator it1 = matchStars.begin() ;  it1!=matchStars.end(); it1++) {
     442    if (!(*it1).ok) continue;
     443    ndata++;
     444    azi[ndata]  = (*it1).azGSC * 3.1415926/180;
     445    elv0[ndata] = (*it1).elvGSC - (*it1).nDiode*1.41/45.;
     446    sig[ndata]  = 0.01;
     447  }
     448   ia[1]   =  ia[2] = 1;
     449   ia[3] = 0;
     450   aa[3] = GondolaGeom::elevSST0;// do not fit elv0
     451
     452  if (ndata<5) return;
     453  float chi2;
     454  try {
     455    lfit(azi, elv0,  sig, ndata, aa, ia, 3, cov, &chi2, sinfunc);
     456  } catch(string s) {
     457    return;
     458  }
     459
     460  double c = aa[1];
     461  double s = aa[2];
     462 
     463  // Get rid of bad fits. The cuts are rather ad hoc
     464 
     465  //if (aa[3] < 39.64 || aa[3] > 39.68) return;
     466  if (chi2/ndata > 4) return;
     467  if (cov[1][1] > 0.0001) return;
     468  if (cov[2][2] > 0.0001) return;
     469
     470  double ampl = sqrt(c*c+s*s);
     471  double phase = atan2(c,s)/(3.1415926/180);
     472 
     473#ifdef STARDUMP
     474  pendstream << snmid << " " << ampl << " " << phase << " " << ndata << " " << chi2/ndata << 
     475       " "  << cov[1][1] << " " << cov[2][2] << '\n';
     476#endif
     477 
     478  pendulInfo info;
     479  info.SN          =   snmid;
     480  info.azPendul    =   phase;
     481  info.angPendul   =   ampl;
     482  pendulInfos[info.SN] = info;
     483
     484  /*
     485  double snum = (matchStars.front().SN + matchStars.back().SN)/2-sn0;
     486  if (snmin > snum || snmax < snum) return;
     487  double elsst = polval(snum, ae, 3);
     488  double azsst = polval(snum, aa, 3);
     489 
     490  if (azsst > 360) azsst -= 360;
     491  if (azsst < 0  ) azsst += 360;
     492  */
     493 
     494//  for (set<TOI>::iterator i = producedTOIs.begin(); i!=producedTOIs.end(); i++) {
     495//    if ((*i).name == "azimuthSST")  computedValue((*i), snum+sn0, azsst);
     496//    if ((*i).name == "elvSST")      computedValue((*i), snum+sn0, elsst);
     497//  }
     498 
     499  free_vector(sn, 1, matchStars.size());
     500  free_vector(elv0, 1, matchStars.size());
     501  free_vector(azi, 1, matchStars.size());
     502  free_vector(sig, 1, matchStars.size());
     503  free_vector(ae, 1, 3);
     504  free_vector(aa, 1, 3);
     505  free_ivector(ia, 1, matchStars.size());
     506  free_matrix(cov, 1, 3, 1, 3);
     507}
     508
     509int StarMatcher::getPendulInfo(double sampleNum, pendulInfo& info) {
     510  map<double, pendulInfo>::iterator i = pendulInfos.lower_bound(sampleNum);
     511  if (i == pendulInfos.begin() && (*i).second.SN >= sampleNum) return -1;
     512  if (i == pendulInfos.end()   && (*i).second.SN <= sampleNum) return -1;
     513 
     514  if ((*i).second.SN > sampleNum) i--;
     515  pendulInfo inf1 = (*i).second;
     516  i++;
     517  pendulInfo inf2 = (*i).second;
     518 
     519  info.SN = sampleNum;
     520  if (inf2.azPendul - inf1.azPendul > 180)  inf2.azPendul -= 360;
     521  if (inf2.azPendul - inf1.azPendul < -180) inf2.azPendul += 360;
     522  info.azPendul  = inf1.azPendul  + (inf2.azPendul  -  inf1.azPendul) * (sampleNum - inf1.SN) / (inf2.SN - inf1.SN);
     523  if (info.azPendul<0)    info.azPendul += 360;
     524  if (info.azPendul>360)  info.azPendul += 360;
     525  info.angPendul = inf1.angPendul + (inf2.angPendul - inf1.angPendul) * (sampleNum - inf1.SN) / (inf2.SN - inf1.SN);
     526  return 0;
    49527}
    50528
    51529
    52530double StarMatcher::getValue(long sampleNum, TOI const& toi) {
    53   if (lastStars.find(sampleNum) == lastStars.end()) return 0;
    54  
    55   SSTEtoile lastStar = lastStars[sampleNum];
    56 
    57 
    58   map<TOI, TOIProducer*> m = neededTOIs[toi];
    59   double lat, lon, ts, alpha, delta, az;
    60  
    61   for (map<TOI, TOIProducer*>::iterator i = m.begin(); i != m.end(); i++) {
    62     TOI inToi = (*i).first;
    63     TOIProducer* prod =  (*i).second;
    64     if (inToi.name == "latitude")  lat = prod->getValue(sampleNum, inToi);
    65     if (inToi.name == "longitude") lon = prod->getValue(sampleNum, inToi);
    66     if (inToi.name == "tsid")       ts = prod->getValue(sampleNum, inToi);
    67     if (inToi.name == "alphaFPC")  alpha = prod->getValue(sampleNum, inToi);
    68     if (inToi.name == "deltaFPC")  delta = prod->getValue(sampleNum, inToi);
    69     if (inToi.name == "azimuthFPC")   az = prod->getValue(sampleNum, inToi);
    70   }
    71    
    72   // find all stars +- 2 deg boresight
    73   double dist = 2;
    74   double dmin = delta - dist; if (dmin<-90) dmin=-90;
    75   double dmax = delta + dist; if (dmax> 90) dmax= 90;
    76   double amin = alpha - dist / cos(delta * 3.1415926/180) / 15.;
    77   if (amin<0) amin += 24;
    78   double amax = alpha + dist / cos(delta * 3.1415926/180) / 15.;
    79   if (amax>24) amax -= 24;
    80  
    81   int a,b,c;
    82   a=0; c=nstars-1;
    83   while (a+1<c) {
    84     b = (a+c)/2;
    85     if (stars[b].dec < dmin) a=b; else c=b;
    86   }
    87   int imin = a;
    88   a=0; c=nstars;
    89   while (a+1<c) {
    90     b = (a+c)/2;
    91     if (stars[b].dec < dmax) a=b; else c=b;
    92   }
    93   int imax = c;
    94 
    95   for (int i=imin; i<=imax; i++) {
    96     if (stars[i].ra >= amin && stars[i].ra <= amax) {
    97       double ha = (ts/3600. - stars[i].ra) * 15. * 3.1415926/180.;
    98       double elv, azim;
    99       hadec_aa(lat * 3.1415926/180., ha, stars[i].dec * 3.1415926/180.,
    100                &elv, &azim);
    101       elv  *= 180/3.1415926;
    102       azim *= 180/3.1415926;
    103       if (azim<0) azim += 360;
    104      
    105       starstream << sampleNum << " " <<
    106         /*lastStar.TEchan << " " <<*/ lastStar.NoDiode << " " <<
    107         alpha << " " << delta << " " <<
    108         az << " " <<
    109         stars[i].ra << " " << stars[i].dec << " " <<
    110         elv << " " << azim << "\n";
    111     }
    112   }
    113  
    114   return 1;
    115 }
     531  processStars();
     532 
     533  // 1. Interpoler la valeur de pendulation
     534  // 2. Interpoler la position en azimuth avec les etoiles encadrant
     535 
     536  pendulInfo pendul;
     537  int rc = getPendulInfo(sampleNum, pendul);
     538  if (rc) return -99999;
     539  if (toi.name == azimuthPendul) return pendul.azPendul;
     540  if (toi.name == anglePendul)   return pendul.angPendul;
     541 
     542  // find nearest matched star
     543  map<double, posInfo>::iterator i = posInfos.lower_bound(sampleNum);
     544  if (i == posInfos.begin() && (*i).second.SN >= sampleNum) return -1;
     545  if (i == posInfos.end()   && (*i).second.SN <= sampleNum) return -1;
     546  if ((*i).second.SN > sampleNum) i--;
     547 
     548  GondolaGeom geom;
     549  geom.setEarthPos((*i).second.lon,(*i).second.lat);
     550  geom.setTSid((*i).second.ts);
     551 
     552  for (map<double, posInfo>::iterator it=i; it != posInfos.end(); it++) {
     553    posInfo s = (*it).second;
     554    double delsn = s.SN - sampleNum;
     555    if (delsn * archParam.acq.perEch > 1) break;
     556    geom.addStar(delsn, s.azStar, s.elvStar, s.diodStar);
     557  }
     558 
     559  if (i != posInfos.begin()) i--;
     560  for (map<double, posInfo>::iterator it=i; it != posInfos.begin(); it--) {
     561    posInfo s = (*it).second;
     562    double delsn = s.SN - sampleNum;
     563    if (-delsn * archParam.acq.perEch > 1) break;
     564    geom.addStar(delsn, s.azStar, s.elvStar, s.diodStar);
     565  }
     566
     567  geom.solveStars();
     568 
     569  if (toi.name == azimuthAxis) return geom.getAzimutAxis();
     570  if (toi.name == elvAxis)     return geom.getElvAxis();
     571  if (toi.name == alphaAxis)   return geom.getAlphaAxis();
     572  if (toi.name == deltaAxis)   return geom.getDeltaAxis();
     573 
     574  if (toi.name == azimuthSST)  return geom.getAzimutSST();
     575  if (toi.name == elvSST)      return geom.getElvSST();
     576  if (toi.name == alphaSST)    return geom.getAlphaSST();
     577  if (toi.name == deltaSST)    return geom.getDeltaSST();
     578 
     579  if (toi.name == azimuthFPC)  return geom.getAzimutCenter();
     580  if (toi.name == elvFPC)      return geom.getElvCenter();
     581  if (toi.name == alphaFPC)    return geom.getAlphaCenter();
     582  if (toi.name == deltaFPC)    return geom.getDeltaCenter();
     583
     584  if (toi.name == azimuthBolo) return geom.getAzimutBolo(toi.index);
     585  if (toi.name == elvBolo)     return geom.getElvBolo(toi.index);
     586  if (toi.name == alphaBolo)   return geom.getAlphaBolo(toi.index);
     587  if (toi.name == deltaBolo)   return geom.getDeltaBolo(toi.index);
     588
     589  return -99999;
     590}
     591
     592bool StarMatcher::canGetValue(long sampleNum, TOI const& /*toi*/) {
     593  processStars();
     594
     595  map<double, pendulInfo>::iterator i = pendulInfos.begin();
     596  if (i == pendulInfos.end()) return false;
     597  if (sampleNum < (*i).second.SN) return false;
     598  i = pendulInfos.end(); i--;
     599  if (sampleNum > (*i).second.SN) return false;
     600 
     601  return true;
     602}
     603
     604bool StarMatcher::canGetValueLater(long sampleNum, TOI const& /*toi*/) {
     605  processStars();
     606 
     607  map<double, pendulInfo>::iterator i = pendulInfos.end();
     608  if (i == pendulInfos.begin()) return true;
     609  i--;
     610  return (sampleNum > (*i).second.SN);
     611}
     612
    116613
    117614
    … …  
    121618  t.insert(TOI("longitude", TOI::unspec, "interp"));
    122619  t.insert(TOI("tsid",      TOI::unspec));
    123   t.insert(TOI("alphaFPC",  TOI::unspec, "galcross0"));
    124   t.insert(TOI("deltaFPC",  TOI::unspec, "galcross0"));
    125   t.insert(TOI("azimuthFPC",  TOI::unspec, "galcross0"));
     620  t.insert(TOI("alphaSST",  TOI::unspec, "galcross0"));
     621  t.insert(TOI("deltaSST",  TOI::unspec, "galcross0"));
     622  t.insert(TOI("azimuthSST",TOI::unspec, "galcross0"));
     623  t.insert(TOI("elvSST",    TOI::unspec, "galcross0"));
     624  t.insert(TOI("rotSpeed",  TOI::unspec, "galcross0"));
    126625  return t;
    127626}
    128627
    129628void StarMatcher::propagateLowBound(TOI const& toi, long sampleNum) {
    130   for (map<long,SSTEtoile>::iterator i = lastStars.begin(); i != lastStars.end(); i++)
    131     if ((*i).first < sampleNum) lastStars.erase(i);
    132 }
    133 
    134 
    135 
    136 
    137 
     629  map<double, posInfo>::iterator i = posInfos.begin();
     630  while (i != posInfos.end() && (*i).first < sampleNum) i++;
     631  if (i != posInfos.begin()) {
     632    i--;
     633    posInfos.erase(posInfos.begin(), i);
     634  }
     635
     636  map<double, pendulInfo>::iterator j = pendulInfos.begin();
     637  while (j != pendulInfos.end() && (*j).first < sampleNum) j++;
     638  if (j != pendulInfos.begin()) {
     639    j--;
     640    pendulInfos.erase(pendulInfos.begin(), j);
     641  }
     642
     643  TOIDerivProducer::propagateLowBound(toi, sampleNum);
     644}
     645
     646
     647// 1. processStars seulement quand au moins 10 etoiles nouvelles
     648// 2. Nettoyer avec fit parabolique sur les 5 dernieres seconde de donnees
     649// 3. Garder periodeRotation secondes de donnees nettoyees
     650// 4. TF ordre 0 sur ces donnees, amplitude et phase -> theta et phi pendulation,
     651//     elevationSST = elv-theta Sin[azimut-phi]
     652//     azimutSST    = azimut+theta Cos[azimut-phi] Tan[elv] (+ OFFSET GALCROSS)
     653//   le signal le plus propre est l'elevation -> fit dessus, puis
     654//   correction azimut SST a partir seconde equation, sans utiliser azimut galcross
     655
     656

  • trunk/Poubelle/archTOI.old/starmatcher.h

    r534 r555  
     1// starmatcher.h
     2// Eric Aubourg         CEA/DAPNIA/SPP   novembre 1999
     3
    14#ifndef STARMATCHER_H
    25#define STARMATCHER_H
    … …  
    47
    58#include "sststarfinder.h"
    6 #include "toipullproducer.h"
     9#include "toiderivproducer.h"
    710
    8 class StarMatcher: public SSTStarProcessor, public TOIPullProducer {
     11class StarMatcher: public SSTStarProcessor, public TOIDerivProducer {
    912public:
    1013  StarMatcher();
    … …  
    1417  virtual string getName();
    1518  virtual double getValue(long sampleNum, TOI const& toi);
     19  virtual bool   canGetValue(long sampleNum, TOI const& toi);
     20  virtual bool   canGetValueLater(long sampleNum, TOI const& toi);
    1621
    1722protected:
    1823  virtual set<TOI>     reqTOIFor(TOI const&);
    19   void  propagateLowBound(TOI const& toi, long sampleNum);
    20 #pragma options align=mac68k
     24  void propagateLowBound(TOI const& toi, long sampleNum);
     25  void processStars();
     26
    2127  struct gscStar {
    2228    float ra;
    … …  
    2430    short mag; // mag * 100
    2531  };
    26 #pragma options align=reset
     32
    2733 
    2834  gscStar* stars;
    2935  long nstars;
    3036 
    31   map<long,SSTEtoile> lastStars;
     37  deque<SSTEtoile> lastStars;
     38 
     39  struct matchStar {
     40    double SN;
     41    double raGSC;
     42    double decGSC;
     43    double azGSC;
     44    double elvGSC;
     45    int    nDiode;
     46    bool   ok;
     47    long   seq;
     48    double lon;
     49    double lat;
     50    double ts;
     51  };
     52 
     53  long lastSeq;
     54 
     55  deque<matchStar> matchStars;
     56  deque<matchStar> cleanStars;
     57 
     58  struct pendulInfo {
     59    double SN;
     60    double azPendul;
     61    double angPendul;
     62  };
     63 
     64  struct posInfo {
     65    double SN;
     66    double azStar;
     67    double elvStar;
     68    double diodStar;
     69    double lon;
     70    double lat;
     71    double ts;
     72  };
     73 
     74  map<double, posInfo> posInfos; // sampleNum -> info
     75  map<double, pendulInfo> pendulInfos; // sampleNum -> info
     76 
     77  int getPendulInfo(double sampleNum, pendulInfo& info); // interpolate
     78 
    3279};
    3380

  • trunk/Poubelle/archTOI.old/templocator.cc

    r534 r555  
    11#include "templocator.h"
     2#include "gondolageom.h"
    23#include <math.h>
    34extern "C" {
    … …  
    161162}
    162163
    163 double TempLocator::getAzimutBolo(int sampleNum, int ibolo) {
    164   double elv, az;
    165   getAltAzBolo(sampleNum, ibolo, elv, az);
     164double TempLocator::getAzimutSST(int sampleNum) {
     165  findGeomFromGC(sampleNum);
     166  double az = azimBolGC;
     167  if (az>360) az -= 360;
     168  if (az<0)   az += 360;
    166169  return az;
    167170}
    168171
    169 double TempLocator::getElvBolo(int sampleNum, int ibolo) {
    170   double elv, az;
    171   getAltAzBolo(sampleNum, ibolo, elv, az);
    172   return elv;
    173 }
    174 
    175 double TempLocator::getAlphaBolo(int sampleNum, int ibolo) {
    176   double elv, az;
    177   getAltAzBolo(sampleNum, ibolo, elv, az);
     172double TempLocator::getElvSST(int /*sampleNum*/) {
     173  return GondolaGeom::elevSST0;
     174}
     175
     176double TempLocator::getAlphaSST(int sampleNum) {
     177  double elv = getElvSST(sampleNum);
     178  double az  = getAzimutSST(sampleNum);
    178179  double ha;
    179180  double ra,dec;
    … …  
    182183  while (ra>24) ra -= 24;
    183184  while (ra<0) ra += 24;
     185  //dec = dec * 180. / M_PI;
     186  return ra;
     187}
     188
     189double TempLocator::getDeltaSST(int sampleNum) {
     190  double elv = getElvSST(sampleNum);
     191  double az  = getAzimutSST(sampleNum);
     192  double ha;
     193  //double ra;
     194  double dec;
     195  aa_hadec (lat * M_PI/180, elv * M_PI/180, az * M_PI/180, &ha, &dec);
     196  //ra = - (ha * 180. / M_PI / 15) + (ts/3600.);
     197  //while (ra>24) ra -= 24;
     198  //while (ra<0) ra += 24;
    184199  dec = dec * 180. / M_PI;
     200  return dec;
     201}
     202
     203
     204
     205double TempLocator::getAzimutBolo(int sampleNum, int ibolo) {
     206  double elv, az;
     207  getAltAzBolo(sampleNum, ibolo, elv, az);
     208  return az;
     209}
     210
     211double TempLocator::getElvBolo(int sampleNum, int ibolo) {
     212  double elv, az;
     213  getAltAzBolo(sampleNum, ibolo, elv, az);
     214  return elv;
     215}
     216
     217double TempLocator::getAlphaBolo(int sampleNum, int ibolo) {
     218  double elv, az;
     219  getAltAzBolo(sampleNum, ibolo, elv, az);
     220  double ha;
     221  double ra,dec;
     222  aa_hadec (lat * M_PI/180, elv * M_PI/180, az * M_PI/180, &ha, &dec);
     223  ra = - (ha * 180. / M_PI / 15) + (ts/3600.);
     224  while (ra>24) ra -= 24;
     225  while (ra<0) ra += 24;
     226  //dec = dec * 180. / M_PI;
    185227  return ra;
    186228}
    … …  
    190232  getAltAzBolo(sampleNum, ibolo, elv, az);
    191233  double ha;
    192   double ra,dec;
    193   aa_hadec (lat * M_PI/180, elv * M_PI/180, az * M_PI/180, &ha, &dec);
    194   ra = - (ha * 180. / M_PI / 15) + (ts/3600.);
     234  //double ra;
     235  double dec;
     236  aa_hadec (lat * M_PI/180, elv * M_PI/180, az * M_PI/180, &ha, &dec);
     237  //ra = - (ha * 180. / M_PI / 15) + (ts/3600.);
     238  //while (ra>24) ra -= 24;
     239  //while (ra<0) ra += 24;
    195240  dec = dec * 180. / M_PI;
    196241  return dec;

  • trunk/Poubelle/archTOI.old/templocator.h

    r534 r555  
    2323  double getAlphaCenter(int sampleNum);
    2424  double getDeltaCenter(int sampleNum);
     25  double getAzimutSST(int sampleNum); // diode 0
     26  double getElvSST(int sampleNum);
     27  double getAlphaSST(int sampleNum);
     28  double getDeltaSST(int sampleNum);
    2529  double getRotSpeed(int sampleNum); // deg/sec
    2630  int    getCrossSamples(int sampleNum, int& SN1, int& SN2);

  • trunk/Poubelle/archTOI.old/toi.cc

    r534 r555  
    6060}
    6161
    62 string TOI::fullName() {
     62string TOI::fullName() const {
    6363 string s = name;
    6464 if (index == all) { s += " all";}
    … …  
    6767   s += string(" ")+idx;
    6868 }
    69  for (set<string>::iterator i = options.begin(); i!=options.end(); i++) {
     69 for (set<string>::const_iterator i = options.begin(); i!=options.end(); i++) {
    7070   s += " " + *i;
    7171 }

  • trunk/Poubelle/archTOI.old/toi.h

    r534 r555  
    3030  void findref();
    3131 
    32   string fullName();
     32  string fullName() const;
    3333};
    3434

  • trunk/Poubelle/archTOI.old/toiboloproducer.cc

    r534 r555  
    4242void TOIBoloProducer::propagateLowBound(TOI const& toi, long sampleNum) {
    4343  CHKPROD
    44   map<TOI, TOIProducer*> need = neededTOIs[toi];
     44  map<TOI, TOIProducer*> & need = neededTOIs[toi];
    4545  int hlen = 1;
    4646  for (map<TOI, TOIProducer*>::iterator i = need.begin(); i != need.end(); i++) {

  • trunk/Poubelle/archTOI.old/toiderivproducer.cc

    r534 r555  
    66#include "archexc.h"
    77#include "requesthandler.h"
     8#include "subsets.h"
    89
    910set<TOI> TOIDerivProducer::reqTOIFor(TOI const&) {
    … …  
    6162  }
    6263 
     64  set<string> unusedextra = extraopts;
     65 
    6366  // 4. Find producers, distribute remaining options
    6467  map<TOI, TOIProducer*> fullInputTOI;
    … …  
    6669  for (set<TOI>::iterator i = inputTOIs.begin(); i != inputTOIs.end(); i++) {
    6770    TOI inTOI = *i;
    68     TOIProducer* prod = TOIManager::findTOIProducer(inTOI);
     71    //TOIProducer* prod0 = TOIManager::findTOIProducer(inTOI);
     72    // Let's see if we get others with options, we have to do this because
     73    // of required options, and we have to iterate over all subsets of the options.
     74    // This algorithm is NOT correct, but should work for realistic cases...
     75    int nmaxopt=-1;
     76    TOIProducer* prod=NULL;
     77    set<set<string> > optsets = getSetOfSubsets(extraopts);
     78    for (set<set<string> >::iterator j = optsets.begin(); j != optsets.end(); j++) {
     79      TOI toi2 = inTOI;
     80      toi2.options.insert((*j).begin(), (*j).end());
     81      TOIProducer* prod0 = TOIManager::findTOIProducer(toi2);
     82      if (prod0 == NULL) continue;
     83      if ((long) toi2.options.size() > nmaxopt) {
     84        prod = prod0;
     85        nmaxopt = toi2.options.size();
     86      }
     87    }
    6988    if (!prod) return false;
    7089    if (!extraopts.empty()) {
    … …  
    7392        if (extraopts.find(*j) != extraopts.end()) {
    7493          inTOI.options.insert(*j);
    75           extraopts.erase(*j);
     94          unusedextra.erase(*j);
    7695        }
    7796      }
    … …  
    8099  }
    81100 
    82   if (!extraopts.empty()) return false;
     101  if (!unusedextra.empty()) return false;
    83102 
    84103  neededTOIs[toi] = fullInputTOI;
    … …  
    91110 // if (!canProduce(toi)) throw ArchExc("cannot produce " + toi.name);
    92111  set<string> s = getProperAvailOptions(toi);
    93   map<TOI, TOIProducer*> need = neededTOIs[toi];
     112  map<TOI, TOIProducer*> & need = neededTOIs[toi];
    94113  for (map<TOI, TOIProducer*>::iterator i = need.begin(); i != need.end(); i++) {
    95114    set<string> s1 = (*i).second->getAvailOptions((*i).first);
    … …  
    101120void TOIDerivProducer::addTOI(TOI& toi, TOIAbsorber* client) {
    102121  TOIProducer::addTOI(toi, client);
    103   map<TOI, TOIProducer*> m = neededTOIs[toi];
     122  map<TOI, TOIProducer*> & m = neededTOIs[toi];
    104123  for (map<TOI, TOIProducer*>::iterator i = m.begin(); i != m.end(); i++) {
    105124    TOI toi2 = (*i).first;
    … …  
    117136void TOIDerivProducer::propagateLowBound(TOI const& toi, long sampleNum) {
    118137  CHKPROD
    119   map<TOI, TOIProducer*> need = neededTOIs[toi];
     138  map<TOI, TOIProducer*> & need = neededTOIs[toi];
    120139  for (map<TOI, TOIProducer*>::iterator i = need.begin(); i != need.end(); i++) {
    121140    (*i).second->wontNeedEarlier((*i).first, this, sampleNum);

  • trunk/Poubelle/archTOI.old/toiinterpolator.cc

    r350 r555  
    11// toiinterpolator.cc
    2 // Eric Aubourg         CEA/DAPNIA/SPP   juillet 1999
     2// Eric Aubourg         CEA/DAPNIA/SPP   octobre 1999
    33
    44
    55#include "toiinterpolator.h"
     6#include "toimanager.h"
     7#include "archexc.h"
    68
    7 void TOIInterpolator::enterValue(double value, long sample)
    8 {
    9   // $CHECK$ sample > back.sample
    10   if (!hist.empty() && sample == hist.back().t) return; // pour eviter /0
    11   hist.push_back(pair(value, sample));
     9TOIInterpolator::TOIInterpolator() {
    1210}
    1311
    14 bool TOIInterpolator::needMoreDataFor(long sample)
    15 {
    16   if (hist.empty()) return true;
    17   int xx = hist.back().t;
    18   return (hist.empty() || sample > xx);
     12
     13string TOIInterpolator::getName() {
     14  return("TOIInterpolator 1.0");
    1915}
    2016
    21 bool TOIInterpolator::canGet(long sample)
    22 {
    23   return (!hist.empty() && (sample >= hist.front().t && sample <= hist.back().t));
     17bool TOIInterpolator::canProduce(TOI const& toi) {
     18  // 1. Already in cache ?
     19  map<TOI, map<TOI, TOIProducer*> >::const_iterator j = neededTOIs.find(toi);
     20  if (j != neededTOIs.end()) return true;
     21
     22  // 2. Should have interp
     23  if (toi.options.find("interp") == toi.options.end()) return false;
     24 
     25  // 3. Can get non interp
     26  TOI toi2 = toi;
     27  toi2.options.erase("interp");
     28  TOIProducer* prod = TOIManager::findTOIProducer(toi2);
     29  set<string> opts = prod->getAvailOptions(toi2);
     30  if (opts.find("interp") != opts.end()) return false; // already handled
     31  if (opts.find("repet") != opts.end()) return false;  // not compatible !
     32 
     33  map<TOI, TOIProducer*> fullInputTOI;
     34  fullInputTOI[toi2] = prod;
     35  neededTOIs[toi] = fullInputTOI;
     36  return true;
    2437}
    2538
    26 double TOIInterpolator::getEValue(long sample) // can forget before sample
    27 {
    28   if (!canGet(sample)) return -9.e99;   // user should have checked before
    29   deque<pair>::iterator i = hist.begin();
    30   while ((*i).t < sample) i++;
    31   if ((*i).t > sample) i--;
    32   double value = (*i).val;
    33   if (i != hist.begin()) {
    34     i--;
    35     hist.erase(hist.begin(), i);
    36   }
    37   return value;
     39set<TOI> TOIInterpolator::reqTOIFor(TOI const& toi) {
     40  set<TOI> x;
     41  if (!canProduce(toi)) return x;
     42  x.insert((*neededTOIs[toi].begin()).first);
     43  return x;
    3844}
    3945
    40 double TOIInterpolator::getIValue(long sample) // can forget before sample
    41 {
    42   if (!canGet(sample)) return -9.e99;   // user should have checked before
     46bool TOIInterpolator::canGetValue(long sampleNum, TOI const& toi) {
     47  map<TOI, TOIProducer*> & inp = neededTOIs[toi];
     48  TOIProducer* prod  = (*inp.begin()).second;
     49  TOI const&   inTOI = (*inp.begin()).first;
    4350 
    44   deque<pair>::iterator i = hist.begin();
    45   while ((*i).t < sample) i++;
    46   double value;
    47   if ((*i).t == sample) {
    48     value = (*i).val;
    49   } else {   
    50     long samp2 = (*i).t;
    51     double val2 = (*i).val;
    52     i--;
    53     long samp1 = (*i).t;
    54     double val1 = (*i).val;
    55     value = val1 + (val2-val1)*(sample-samp1)/(samp2-samp1);
    56   }
    57   if (i != hist.begin()) {
    58     i--;
    59     hist.erase(hist.begin(), i);
    60   }
    61   return value;
    62 }
    63 
    64 bool TOIInterpolator::isNewValue(long sample)
    65 {
    66   deque<pair>::iterator i = hist.begin();
    67   while (i!=hist.end()) {
    68     if ((*i).t == sample) return true;
    69     if ((*i).t > sample) return false;
    70     i++;
    71   }
     51  if (prod->canGetValue(sampleNum, inTOI)) return true; // direct
     52  if (prod->canGetPrevValue(sampleNum, inTOI) &&
     53      prod->canGetNextValue(sampleNum, inTOI)) return true;
    7254  return false;
    7355}
    7456
    75 
    76 long TOIInterpolator::nextSample(long sampleMin)
    77 {
    78   for (deque<pair>::iterator i = hist.begin(); i<hist.end(); i++)
    79     if ((*i).t >= sampleMin) return (*i).t;
    80   return -1;
     57bool TOIInterpolator::canGetValueLater(long sampleNum, TOI const& toi) {
     58  map<TOI, TOIProducer*> & inp = neededTOIs[toi];
     59  TOIProducer* prod  = (*inp.begin()).second;
     60  TOI const&   inTOI = (*inp.begin()).first;
     61 
     62  if (prod->canGetValue(sampleNum, inTOI)) return false; // direct
     63  if (prod->canGetValueLater(sampleNum, inTOI)) return true; // direct
     64  if (!prod->canGetPrevValue(sampleNum, inTOI)) return false; // no hope
     65  if (prod->canGetNextValue(sampleNum, inTOI)) return false; // can get now
     66  return true;
    8167}
    8268
     69double TOIInterpolator::getValue(long sampleNum, TOI const& toi) {
     70  map<TOI, TOIProducer*> & inp = neededTOIs[toi];
     71  TOIProducer* prod  = (*inp.begin()).second;
     72  TOI const&   inTOI = (*inp.begin()).first;
     73
     74  if (prod->canGetValue(sampleNum, inTOI))
     75     return prod->getValue(sampleNum, inTOI); // direct
     76     
     77  long sn1 = sampleNum;
     78  double v1 = prod->getPrevValue(sn1, inTOI);
     79 
     80  long sn2 = sampleNum;
     81  double v2 = prod->getNextValue(sn2, inTOI);
     82
     83  if (sn1 == sn2) throw ArchExc("Inconsistence in interp, sn1==sn2");
     84
     85  return v1 + (v2-v1) * (sampleNum-sn1) / (sn2-sn1);
     86}
     87
     88
     89void TOIInterpolator::propagateLowBound(TOI const& toi, long sampleNum) {
     90  CHKPROD
     91  map<TOI, TOIProducer*> & inp = neededTOIs[toi];
     92  TOIProducer* prod  = (*inp.begin()).second;
     93  TOI const&   inTOI = (*inp.begin()).first;
     94  if (prod->canGetPrevValue(sampleNum,toi)) {
     95    prod->getPrevValue(sampleNum, toi);
     96    sampleNum--;
     97    prod->wontNeedEarlier(toi, this, sampleNum);
     98  }
     99}
     100

  • trunk/Poubelle/archTOI.old/toiinterpolator.h

    r350 r555  
    11// toiinterpolator.h
    2 // Eric Aubourg         CEA/DAPNIA/SPP   juillet 1999
     2// Eric Aubourg         CEA/DAPNIA/SPP   octobre 1999
    33
    44
    … …  
    66#define TOIINTERPOLATOR_H
    77
    8 #include <deque>
    9 
    10 using namespace std;
     8#include "toipullproducer.h"
    119
    1210
    13 class TOIInterpolator {
     11class TOIInterpolator : public TOIPullProducer {
    1412public:
    15 
    16   void enterValue(double value, long sample);
    17   bool   needMoreDataFor(long sample);   // need more data ? (for interp)
    18   bool   canGet(long sample);        // enough info ?
    19   double getIValue(long sample);     // can forget before sample
    20   double getEValue(long sample);     // can forget before sample
    21   bool   isNewValue(long sample);   
    22   long   nextSample(long sampleMin); // the next true sample available
     13  TOIInterpolator();
     14  virtual string getName();
     15  virtual double getValue(long sampleNum, TOI const& toi);
     16  virtual bool   canGetValue(long sampleNum, TOI const& toi);
     17  virtual bool   canGetValueLater(long sampleNum, TOI const& toi);
     18  virtual bool   canProduce(TOI const&);
    2319
    2420protected:
    25   struct pair {
    26     pair(double v, long tt) : val(v), t(tt) {}
    27     pair() {}
    28     double val;
    29     long t;
    30   };
    31  
    32   deque<pair> hist;
     21  virtual set<TOI>     reqTOIFor(TOI const&);
     22  void                 propagateLowBound(TOI const& toi, long sampleNum);
     23
    3324};
    3425

  • trunk/Poubelle/archTOI.old/toiiter.cc

    r534 r555  
    272272  if (curSample >= sEnd || archParam.acq.SN2MJD(curSample) >= mjdEnd) return false;
    273273 
    274   for (int i=0; i<request.size(); i++)
    275     request[i].second->wontNeedEarlier(request[i].first, this, curSample);
     274  static long lastClean = 0;
     275  if (curSample - lastClean > 100) {
     276    for (int i=0; i<request.size(); i++)
     277      request[i].second->wontNeedEarlier(request[i].first, this, curSample);
     278    lastClean = curSample;
     279  }
    276280
    277281  bool endFound = false;

  • trunk/Poubelle/archTOI.old/toimanager.cc

    r534 r555  
    7777
    7878#include "timetoiproducer.h"
     79#include "toiinterpolator.h"
     80#include "toirepeater.h"
     81#include "toiflagger.h"
    7982#include "toiboloproducer.h"
    8083#include "sststarfinder.h"
    … …  
    8285#include "galcrosslocator.h"
    8386#include "starmatcher.h"
     87#include "rotspeed.h"
     88#include "gyrocalibrator.h"
    8489
    8590void TOIManager::registerDefaultProducers() {
    … …  
    9297 
    9398  REGPROD(TimeTOIProducer);
     99  REGPROD(TOIInterpolator);
     100  REGPROD(TOIRepeater);
     101  REGPROD(TOIFlagger);
    94102  REGPROD(TOIBoloProducer);
    95103  REGPROD(SSTStarFinder);
    … …  
    97105  REGPROD(GalCrossLocator);
    98106  REGPROD(StarMatcher);
     107  REGPROD(RotSpeed);
     108  REGPROD(GyroCalibrator);
    99109}
    100110

  • trunk/Poubelle/archTOI.old/toiproducer.cc

    r534 r555  
    5757  // should be nothing left
    5858  if (!extraopts.empty()) return false;
    59  
     59
    6060  return true;
    6161}
    … …  
    129129
    130130class comp2nd { public:bool operator()(pair<long, double> const& a, long b) {return a.first<b;}};
    131 
    132 deque<pair<long, double> >::const_iterator TOIProducer::findHist(TOI const& toi,long sampleNum) {
     131class comp2ndI { public:bool operator()(long a, pair<long, double> const& b) {return a<b.first;}};
     132
     133deque<pair<long, double> >::const_iterator TOIProducer::findHistL(TOI const& toi,long sampleNum) {
    133134  CHKPROD
    134135  deque<pair<long, double> >& h = history[toi.ref];
    … …  
    136137}
    137138
    138 deque<pair<long, double> >::const_iterator TOIProducer::findHist(deque<pair<long, double> >& h,long sampleNum) {
     139deque<pair<long, double> >::const_iterator TOIProducer::findHistL(deque<pair<long, double> >& h,long sampleNum) {
    139140  return lower_bound(h.begin(), h.end(), sampleNum, comp2nd());
    140141}
     142
     143deque<pair<long, double> >::const_iterator TOIProducer::findHistH(TOI const& toi,long sampleNum) {
     144  CHKPROD
     145  deque<pair<long, double> >& h = history[toi.ref];
     146  return upper_bound(h.begin(), h.end(), sampleNum, comp2ndI());
     147}
     148
     149deque<pair<long, double> >::const_iterator TOIProducer::findHistH(deque<pair<long, double> >& h,long sampleNum) {
     150  return upper_bound(h.begin(), h.end(), sampleNum, comp2ndI());
     151}
     152
    141153
    142154bool TOIProducer::canGetValue(long sampleNum, TOI const& toi) {
    … …  
    145157  if (h.empty()) return false;
    146158  if (sampleNum < h.front().first || sampleNum > h.back().first) return false;
    147   deque<pair<long, double> >::const_iterator i = findHist(h, sampleNum);
     159  deque<pair<long, double> >::const_iterator i = findHistL(h, sampleNum);
    148160  if (i == h.end()) return false;
    149161  return ((*i).first == sampleNum);
    … …  
    163175  deque<pair<long, double> >& h = history[toi.ref];
    164176  if (h.empty()) return false;
    165   return (sampleNum > h.back().first && sampleNum <= h.front().first);
     177  return (sampleNum > h.front().first && sampleNum <= h.back().first);
    166178   // Must be inside, otherwise we don't know if it is prev or simply past
    167179}
    … …  
    171183  deque<pair<long, double> >& h = history[toi.ref];
    172184  if (h.empty()) return false;
    173   return (sampleNum >= h.back().first && sampleNum < h.front().first);
     185  return (sampleNum >= h.front().first && sampleNum < h.back().first);
    174186}
    175187
    … …  
    179191//       << h.front().first
    180192//       << " " << h.back().first << endl;
    181   deque<pair<long, double> >::const_iterator i = findHist(h, sampleNum);
     193  deque<pair<long, double> >::const_iterator i = findHistL(h, sampleNum);
    182194  if (i == h.end()) return -999999;
    183195//  cout << "found " << (*i).first << endl;
    … …  
    189201double TOIProducer::getPrevValue(long& sampleNum, TOI const& toi) {
    190202  deque<pair<long, double> >& h = history[toi.ref];
    191   deque<pair<long, double> >::const_iterator i = findHist(h, sampleNum);
     203  deque<pair<long, double> >::const_iterator i = findHistH(h, sampleNum);
    192204  if (i == h.end()) return -999999;
    193   if ((*i).first == sampleNum) {
     205  pair<long, double> tstpair = *i;
     206  while ((*i).first >= sampleNum) {
    194207    if (i == history[toi.ref].begin()) {sampleNum = -99999 ;return -9999;}
    195208    i--;
    196209  }
     210  tstpair = *i;
    197211  sampleNum = (*i).first;
    198212  return (*i).second;
    … …  
    201215double TOIProducer::getNextValue(long& sampleNum, TOI const& toi) {
    202216  deque<pair<long, double> >& h = history[toi.ref];
    203   deque<pair<long, double> >::const_iterator i = findHist(h, sampleNum);
     217  deque<pair<long, double> >::const_iterator i = findHistL(h, sampleNum);
    204218  if (i == h.end()) return -999999;
    205   if ((*i).first == sampleNum) {
     219  pair<long, double> tstpair = *i;
     220  while ((*i).first <= sampleNum) {
    206221    if (i == history[toi.ref].end()) {sampleNum = -99999 ;return -9999;}
    207222    i++;
    208223  }
     224  tstpair=*i;
    209225  sampleNum = (*i).first;
    210226  return (*i).second;
    … …  
    212228
    213229void TOIProducer::computedValue(TOI const& toi, long sampleNum, double value) {
    214   map<int, deque<pair<long, double> > >::iterator i = history.find(toi.ref);
    215   if (i == history.end()) throw ArchExc("computedValue : bad TOI " + toi.name);
    216   deque<pair<long, double> >& h = (*i).second;
     230  //map<int, deque<pair<long, double> > >::iterator i = history.find(toi.ref);
     231  //if (i == history.end()) throw ArchExc("computedValue : bad TOI " + toi.name);
     232  //deque<pair<long, double> >& h = (*i).second;
     233  deque<pair<long, double> >& h = history[toi.ref];
    217234  if (!h.empty() && sampleNum <= h.back().first)
    218235    throw ArchExc("computedValue : sampleNum not in sequence for " + toi.name);
    219236  h.push_back(pair<long,double>(sampleNum,value));
    220  
    221237  for (map<TOIAbsorber*, set<TOI> >::iterator j = clients.begin(); j != clients.end(); j++) {
    222238    set<TOI>& tois = (*j).second;

  • trunk/Poubelle/archTOI.old/toiproducer.h

    r534 r555  
    6060  map<int, deque<pair<long, double> > >  history;
    6161 
    62   deque<pair<long, double> >::const_iterator findHist(TOI const& toi, long sampleNum);
    63   deque<pair<long, double> >::const_iterator findHist(deque<pair<long, double> >& h, long sampleNum);
     62  deque<pair<long, double> >::const_iterator findHistL(TOI const& toi, long sampleNum);
     63  deque<pair<long, double> >::const_iterator findHistL(deque<pair<long, double> >& h, long sampleNum);
     64  deque<pair<long, double> >::const_iterator findHistH(TOI const& toi, long sampleNum);
     65  deque<pair<long, double> >::const_iterator findHistH(deque<pair<long, double> >& h, long sampleNum);
    6466};
    6567

  • trunk/Poubelle/archTOI.old/toipullproducer.cc

    r534 r555  
    1010
    1111bool TOIPullProducer::canGetValue(long sampleNum, TOI const& toi) {
    12   map<TOI, TOIProducer*> m = neededTOIs[toi];
     12  map<TOI, TOIProducer*> & m = neededTOIs[toi];
    1313  for (map<TOI, TOIProducer*>::iterator i = m.begin(); i != m.end(); i++) {
    1414    if (!(*i).second->canGetValue(sampleNum-needBefore, (*i).first) ||
    … …  
    2727
    2828bool TOIPullProducer::canGetValueLater(long sampleNum, TOI const& toi) {
    29   map<TOI, TOIProducer*> m = neededTOIs[toi];
     29  map<TOI, TOIProducer*> & m = neededTOIs[toi];
    3030  for (map<TOI, TOIProducer*>::iterator i = m.begin(); i != m.end(); i++) {
    3131    if (!(*i).second->canGetValueLater(sampleNum+needAfter, (*i).first)) return false;
    … …  
    4646long TOIPullProducer::firstSampleNum(TOI const& toi) {
    4747  long xx = -999999999L;
    48   map<TOI, TOIProducer*> m = neededTOIs[toi];
     48  map<TOI, TOIProducer*> & m = neededTOIs[toi];
    4949  for (map<TOI, TOIProducer*>::iterator i = m.begin(); i != m.end(); i++) {
    5050    long x = (*i).second->firstSampleNum((*i).first) + needBefore;
    … …  
    5656long TOIPullProducer::lastSampleNum(TOI const& toi) {
    5757  long xx = 999999999L;
    58   map<TOI, TOIProducer*> m = neededTOIs[toi];
     58  map<TOI, TOIProducer*> & m = neededTOIs[toi];
    5959  for (map<TOI, TOIProducer*>::iterator i = m.begin(); i != m.end(); i++) {
    6060    long x = (*i).second->lastSampleNum((*i).first) - needAfter;
    … …  
    6666void TOIPullProducer::propagateLowBound(TOI const& toi, long sampleNum) {
    6767  CHKPROD
    68   map<TOI, TOIProducer*> need = neededTOIs[toi];
     68  map<TOI, TOIProducer*> & need = neededTOIs[toi];
    6969  for (map<TOI, TOIProducer*>::iterator i = need.begin(); i != need.end(); i++) {
    7070    (*i).second->wontNeedEarlier((*i).first, this, sampleNum-needBefore);

  • trunk/Poubelle/archTOI.old/tsidproducer.cc

    r534 r555  
    2222
    2323double TSidProducer::getValue(long sampleNum, TOI const& toi) {
    24   map<TOI, TOIProducer*> m = neededTOIs[toi];
     24  map<TOI, TOIProducer*> & m = neededTOIs[toi];
    2525  // Seulement longitude...
    26   TOI longTOI = (*m.begin()).first;
     26  TOI const& longTOI = (*m.begin()).first;
    2727  TOIProducer* longProd =  (*m.begin()).second;
    2828  myTS.setLongitude(longProd->getValue(sampleNum, longTOI));
Note: See TracChangeset for help on using the changeset viewer.