source: Sophya/trunk/Poubelle/archTOI.old/starmatcher.cc@ 555

Last change on this file since 555 was 555, checked in by ansari, 26 years ago

portage cxx en cours
File size: 20.8 KB
Line 
1// starmatcher.cc
2// Eric Aubourg CEA/DAPNIA/SPP novembre 1999
3
4#include "starmatcher.h"
5#include "sststarfinder.h"
6#include "toimanager.h"
7#include "archexc.h"
8#include "archparam.h"
9#include "gondolageom.h"
10
11extern "C" {
12#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;
44}
45
46#include <stdio.h>
47
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
95StarMatcher::StarMatcher() {
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
115 FILE* f;
116
117 f = fopen("gsc7.dat","r");
118 if (!f) throw ArchExc("Error opening gsc7.dat");
119
120 int4 n4;
121 fread(&n4, sizeof(int4), 1 , f);
122
123#ifdef SWAP
124 bswap4(&n4);
125#endif
126 nstars = n4;
127
128 stars = new gscStar[nstars];
129 char* compdata = new char[10*nstars];
130 fread(compdata, 10, nstars, f);
131 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;
160
161 TOIProducer* prod = TOIManager::findTOIProducer(TOI("sstStarCount"));
162 if (!prod) {
163 cerr << "StarMatcher : cannot find producer for sstStarCount" << endl;
164 exit(-1);
165 }
166
167 SSTStarFinder* sprod = dynamic_cast<SSTStarFinder*>(prod);
168 if (!sprod) {
169 cerr << "StarMatcher : producer for sstStarCount is not a SSTStarFinder" << endl;
170 exit(-1);
171 }
172
173 lastSeq = 0;
174
175 sprod->registerProcessor(this);
176
177}
178
179string StarMatcher::getName() {
180 return("StarMatcher 1.0");
181}
182
183#ifdef STARDUMP
184static ofstream starstream("stars.dat");
185static ofstream cstarstream("cstars.dat");
186static ofstream pendstream("pendul.dat");
187#endif
188
189void StarMatcher::dataFeed(SSTEtoile const& 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;
527}
528
529
530double StarMatcher::getValue(long sampleNum, TOI const& toi) {
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
613
614
615set<TOI> StarMatcher::reqTOIFor(TOI const&) {
616 set<TOI> t;
617 t.insert(TOI("latitude", TOI::unspec, "interp"));
618 t.insert(TOI("longitude", TOI::unspec, "interp"));
619 t.insert(TOI("tsid", TOI::unspec));
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"));
625 return t;
626}
627
628void StarMatcher::propagateLowBound(TOI const& toi, long sampleNum) {
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
Note: See TracBrowser for help on using the repository browser.