source: Sophya/trunk/Poubelle/archediab.old/SimulMission/simulmission.c@ 647

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

simulation
File size: 28.6 KB
Line 
1/* Simulation de donnees mission, Eric Aubourg & Jacques Delabrouille, mai 1999 */
2
3//#include <stdio.h>
4//#include <string.h>
5#include <math.h>
6#include <stdlib.h>
7
8#include "manip.h"
9#include "choix_acquisition.h"
10#include "archeops.h"
11#include "choix_param.h"
12#include "structure.h"
13#include "tm.h"
14#include "simulmission.h"
15#include "simulstate.h"
16#include "compress.h"
17
18#define M_PI 3.14159265358979
19
20double GauRnd(double am, double s);
21
22
23#define drand01() ( (double)rand() / RAND_MAX )
24
25
26
27
28
29
30static int numblock = 0;
31
32static int stateChanged = 0;
33
34static long tickStart = -1;
35
36double overSpeed = 1;
37double whiteNoise[6];
38double bolfreq[6];
39double glitchFreq; // Hz
40double glitchMaxAmpl; // W
41double boloTimeCst; // s
42int sst2Bars;
43
44
45param_bolo simParam;
46reglage_bolo simReglage;
47
48// Fichier avec donnees mission
49static short fmisdat = 0;
50
51struct boloobs {
52 double power;
53};
54
55struct posbal {
56 double lon; // degres
57 double lat; // degres
58 double azimut; // degres
59};
60
61struct misdat {
62 double mjd; // JD - 2450000 --- 9 octobre 1995 midi UTC
63 struct posbal pos; // Position du ballon
64 struct boloobs bolo[6]; // Puissance sur chaque bolo
65};
66
67
68struct localgeom {
69 double ra; // Centre du pointage
70 double dec;
71 double dra; // vecteur de deplacement, normalise a une minute d'arc
72 double ddec;
73 double ora; // vecteur orthogonal au deplacement, vers le zenith
74 double odec; // normalise a une minute d'arc
75};
76
77static int szmisdatbuf;
78static struct misdat* misdatbuf[2];
79
80static int curBuf = 0;
81static int iBuf = 0;
82static int needBuf = -1;
83
84static int iRdsample; // celui de lastDat
85static int iGensample; // 0..oversample-1
86static int oversample; // surechantillonnage necessaire...
87
88// Donnees GSC pour simulation senseur stellaire
89
90struct star {
91 float ra;
92 float dec;
93 short mag; // mag * 100
94};
95
96static struct star* stars=0; // trie par dec croissant...
97long nstars = 0;
98
99// Donnees correspondant a un scan
100long databolo[nb_max_bolo][nb_per_block*2];
101int datadiod[nb_per_block*2][48];
102double latitude, longitude, mjd;
103
104static double per_echant(void);
105static void restartTime(void);
106
107#define lastData misdatbuf[curBuf][iBuf]
108#define nextData misdatbuf[curBuf][iBuf+1]
109
110// Valeur precedente du bolometre pour gestion inertie (glitchs...) en muV
111static double boloLast[nb_max_bolo];
112
113// Diverses fonctions utiles...
114void calcLocalGeom(struct misdat* md, struct localgeom* geom);
115
116
117static void restartTime(void) {
118 tickStart = TickCount() - numblock * per_echant() * (double)nb_per_block*2.*60.;
119}
120
121static volatile int hasDebugInfo=0;
122static volatile char debugInfo[80];
123
124static long readblocks = 0;
125
126extern param_bolo parametr; // dans archeops.c
127extern reglage_bolo reglage_standard[8]; // liste bolo dans le programme principal
128
129static double laststarmag[46][10];
130static double laststartime[46][10];
131
132void InitSimulMission(void) {
133 long n,nn;
134 long nSamples;
135
136 for (n=0; n<6; n++) {
137 bolfreq[n] = 100e9;
138 whiteNoise[n] = 2.e-17;
139 boloLast[n] = 0;
140 }
141
142 for (n=0; n<46; n++)
143 for (nn=0; nn<10; nn++) {
144 laststarmag[n][nn] = laststartime[n][nn] = 0;
145 }
146
147 glitchFreq = 0.3;
148 glitchMaxAmpl = 1.e-11;
149 boloTimeCst = 1.e-2;
150
151 // Datacards
152 ReadSimulDC();
153
154 // Donnees mission, initialisation du buffer
155 if (FSOpen("\pmission.dat", 0, &fmisdat)) {
156 printf("***ERREUR*** ouverture mission.dat\n");
157 return;
158 }
159
160 szmisdatbuf = 512;
161 misdatbuf[0] = (struct misdat*) malloc(szmisdatbuf * sizeof(struct misdat));
162 misdatbuf[1] = (struct misdat*) malloc(szmisdatbuf * sizeof(struct misdat));
163 curBuf = 0;
164 iBuf=0;
165 needBuf = -1;
166 iRdsample = 0;
167 iGensample = 0;
168
169 n = sizeof(long);
170 FSRead(fmisdat, &n, (char*) &nSamples);
171 printf("Fichier mission.dat : %d echantillons\n", nSamples);
172 n = szmisdatbuf * sizeof(struct misdat);
173 FSRead(fmisdat, &n, (char*) misdatbuf[0]);
174 n = (szmisdatbuf-1) * sizeof(struct misdat);
175 memcpy(misdatbuf[1], &misdatbuf[0][szmisdatbuf-1], sizeof(struct misdat));
176 FSRead(fmisdat, &n, (char*) misdatbuf[1]);
177 readblocks = 2*szmisdatbuf;
178 // Nb d'echantillons dans un cercle, nominalement : 1 cercle = 30 secondes,
179 // ech a 180 Hz, et en realite 125 ech par cercle dans le fichier...
180 oversample = (int) ((30 * 180.) / 125.);
181 printf("oversample = %d\n",oversample);
182
183 // Catalogue GSC pour senseur stellaire
184 {
185 short ref;
186 long n;
187 //if (FSOpen("\phyp6HADY.dat",0,&ref)) {
188 // printf("***ERREUR*** ouverture hyp6HADY.dat, on va utiliser gsc7\n");
189 if (FSOpen("\pgsc7.dat",0,&ref)) {
190 printf("***ERREUR*** ouverture gsc7.dat\n");
191 return;
192 }
193 //}
194 n = sizeof(long);
195 FSRead(ref, &n, (char*) &nstars);
196 n = nstars*sizeof(struct star);
197 stars = (struct star*) malloc(n);
198 FSRead(ref, &n, (char*) stars);
199 FSClose(ref);
200 printf("GSC m<7 : read %d stars\n", nstars);
201 }
202 hasDebugInfo=0;
203 simReglage = reglage_standard[0];
204 simParam = parametr;
205}
206
207
208void SimulMissionReadLoop(void)
209{
210 long n;
211 short rc;
212 if (hasDebugInfo) {
213 printf("*** Debug **%d**\n",hasDebugInfo);
214 hasDebugInfo=0;
215 }
216 if (needBuf < 0) return;
217 if (fmisdat == 0) return;
218 readblocks += szmisdatbuf;
219
220 printf("Reading simulation data block %d\n",readblocks);
221 memcpy(misdatbuf[needBuf], &misdatbuf[1-needBuf][szmisdatbuf-1], sizeof(struct misdat));
222 n = (szmisdatbuf-1) * sizeof(struct misdat);
223 if ((rc=FSRead(fmisdat, &n, (char*) misdatbuf[needBuf]))!=0) {
224 printf("Erreur lecture %d\n", rc);
225 FSClose(fmisdat);
226 fmisdat = 0;
227 }
228 needBuf = -1;
229}
230
231float isStarInArray(float ra, float dec, float alow, float dlow,
232 float a1, float d1, float a2, float d2) // -1 = no, 0-1 : yes
233{
234 float x,y;
235 ra = ra-alow;
236 if (ra>12) ra -=24;
237 if (ra<-12) ra += 24;
238 x = (15*cos(dec * 3.1415926/180));
239 ra *= x;
240 a1 *= x;
241 a2 *= x;
242 x = ((ra)*a1 + (dec-dlow)*d1)/(a1*a1 + d1*d1);
243 y = ((ra)*a2 + (dec-dlow)*d2)/(a2*a2 + d2*d2);
244 if (x<0 || x>1 || y<0 || y>1) return -1;
245 return y;
246}
247
248// Constantes de temps du senseur stellaire
249static const double t1 = 17e-3;
250static const double t2 = 150e-3;
251static const double t3 = 0;
252static const double dt = 30e-3;
253
254static double sig1(double t);
255static double sig2(double t);
256
257static double sig1(double t)
258{
259 return (t2*(t1-t3)*exp(-t/t1) - t1*(t2-t3)*exp(-t/t2))/(t1-t2)/t1;
260}
261
262static double sig2(double t)
263{
264 double x=0;
265 if (t>=0) x -= sig1(t);
266 if (t>=dt) x += sig1(t-dt);
267 return x;
268}
269
270// Balayage dans la direction dra, ddec, pointage central sur ra, dec
271void SSTSignal(float ra, float dec, float dra, float ddec,
272 float ora, float odec, int* diodes, double secondes)
273{
274 // Valeurs precedentes pour gestion des constantes de temps
275
276 //float ora, odec; // orthogonal to dra, ddec, directed to zenith
277 float dmin, dmax; // min and max value of dec in which to look for stars
278 int imin, imax; // corresponding indices in gsc catalog
279 float dlow, alow; // "bottom left" corner of first diode array
280 float d1, a1, d2, a2; // sides of diode array
281 float dd, dr; // translation from first to second diode array.
282 int a,b,c;
283 // renormalize dra, ddec to have a norm of 1 arc min on sky
284 {
285 double nrm;
286 nrm = dra*15*cos(dec * 3.1415926/180);
287 nrm = sqrt(ddec*ddec + nrm*nrm)*60;
288 ddec /= nrm;
289 dra /= nrm;
290 }
291 // orthogonal, vertical, vector. up-down depends on spin orientation.
292 //ora = ddec / (15*cos(dec * 3.1415926/180));
293 //odec = - dra * (15*cos(dec * 3.1415926/180));
294
295 // size of diode array
296 // along spin, 4*1.88 arc minutes
297 a1 = 4*1.88*dra; d1 = 4*1.88*ddec;
298 // orthogonal, 46 diodes of 1.88 arc minutes each
299 // (mail from Silva Masi, 1.41 deg in elevation between center of pixel 1
300 // and center of pixel 46).
301 a2 = (1.88*46)*ora; d2 = (1.88*46)*odec;
302
303 alow = ra - a1/2 - a2/2;
304 dlow = dec - d1/2 - d2/2;
305
306 // shift with second diode array
307 dd = 10. * ddec - .5 * odec;
308 dr = 10. * dra - .5 * ora;
309
310 // Fin min/max values of delta, and thus imin imax in gsc..
311 dmin = dmax = dlow;
312 if ((dlow+d1) < dmin) dmin = dlow+d1;
313 if ((dlow+d2) < dmin) dmin = dlow+d2;
314 if ((dlow+d1+d2) < dmin) dmin = dlow+d1+d2;
315 if ((dlow+d1) > dmax) dmax = dlow+d1;
316 if ((dlow+d2) > dmax) dmax = dlow+d2;
317 if ((dlow+d1+d2) > dmax) dmax = dlow+d1+d2;
318 if ((dmin+dd) < dmin) dmin = dmin+dd;
319 if ((dmax+dd) > dmax) dmax = dmax+dd;
320
321 a=0; c=nstars-1;
322 while (a+1<c) {
323 b = (a+c)/2;
324 if (stars[b].dec < dmin) a=b; else c=b;
325 }
326 imin = a;
327 a=0; c=nstars;
328 while (a+1<c) {
329 b = (a+c)/2;
330 if (stars[b].dec < dmax) a=b; else c=b;
331 }
332 imax = c;
333
334 //for (c=0; c<46; c++) diodes[c] = 50; // fond = 50...
335
336 for (c=imin; c<=imax; c++) {
337 float f = isStarInArray(stars[c].ra, stars[c].dec,
338 alow, dlow, a1, d1, a2, d2);
339 // Pas += a cause de doublons dans le GSC. On n'a jamais plus d'une etoile...
340 //if (f>=0) diodes[(int)(f*46)] = 1000-stars[c].mag;
341 if (f>=0) {
342 int idiod = (int)(f*46);
343 double flux = 3500*pow(10,-(stars[c].mag/250.));
344 if (secondes - laststartime[idiod][9]>dt/3 || laststarmag[idiod][9] != flux) {
345 // mag=0 -> flux = 3500 (Vega)
346 // decalage...
347 for (a=0; a<9; a++) {
348 laststartime[idiod][a] = laststartime[idiod][a+1];
349 laststarmag[idiod][a] = laststarmag[idiod][a+1];
350 }
351 laststartime[idiod][9] = secondes;
352 laststarmag[idiod][9] = flux;
353 }
354 }
355
356 // Test dans la deuxieme barrette
357 if (sst2Bars) {
358 f = isStarInArray(stars[c].ra, stars[c].dec,
359 alow+dr, dlow+dd, a1, d1, a2, d2);
360 //if (f>=0) diodes[(int)(f*46)] = 1000-stars[c].mag;
361 if (f>=0) {
362 int idiod = (int)(f*46);
363 double flux = 3500*pow(10,-(stars[c].mag/250.));
364 if (secondes - laststartime[idiod][9]>dt/3 || laststarmag[idiod][9] != flux) {
365 for (a=0; a<9; a++) {
366 laststartime[idiod][a] = laststartime[idiod][a+1];
367 laststarmag[idiod][a] = laststarmag[idiod][a+1];
368 }
369 laststartime[idiod][9] = secondes;
370 laststarmag[idiod][9] = flux;
371 }
372 }
373 }
374 }
375 for (c=0; c<46; c++) {
376 double signal = 0;
377 for (a=0; a<10; a++)
378 signal += sig2(secondes - laststartime[c][a]) * laststarmag[c][a];
379 signal += (c*2-46); // offset specifique a chaque diode
380 if (signal > 2047) signal = 2047;
381 if (signal < -2047) signal = -2047;
382 //if (signal < 0) signal += 4096;
383 signal += 2048;
384 diodes[c] = signal;
385 //if ( (secondes - laststartime[c])<dt && c>0) diodes[c-1] = laststarmag[c];
386 //diodes[c] = (secondes - laststartime[c])<dt ? laststarmag[c] : 0;
387 //diodes[c] = laststarmag[c] ;
388 }
389}
390
391static int vitesse[20] = {10, // param
392 0, // journal
393 10, // reglage
394 0, // dilution
395 10, // gps
396 10, // 1 periode
397 0, // synchro sol
398 0, // pointage sol
399 0, // bolo
400 0, // gyro
401 1, // sst
402 0, // 11
403 1, // bolo comp
404 0, // gyro comp
405 0, // sst comp
406 0, // status flash
407 0, // cmd flash
408 0, // data brut
409 0, // 18
410 0}; // 19
411
412static double per_echant(void)
413{
414 double p,f1,f2,f3;
415 int pp;
416 pp=simReglage.horloge.periode;
417 p=pp/5.;
418 f1=1000/p;f2=f1/bitmot;f3=f2*1000./(double)(simReglage.horloge.nb_mesures);
419// printf(" per_echantillon : srhp %f srhn %f f3 %f \n",pp,simReglage.horloge.nb_mesures,f3);
420 return 0.5/f3; // 2 fois la frequence de modulation
421}
422
423int emission_blocks(int type);
424
425extern int mode_transmission_telemesure[nb_modes_telemesure][nb_type_blocks];
426
427
428int emission_blocks(int type)
429{
430int a,t;//,e;
431t = (simReglage.dilu.transmission&0xf);
432//e = (simReglage.dilu.transmission&0x30)>>4;
433t = mode_transmission_telemesure[t][type];
434
435//e = mode_transmission_flash[e][type];
436a=0;
437if( (t >0) && ( (numblock%t)==0) ) a=1;
438//if( (e >0) && ( (numblock%e)==0) ) a+=2;
439return(a);
440}
441
442
443int SimulMissionBloc(tmtc* tt)
444{
445 double t;
446 static int iLast = -99;
447
448 if (tickStart<0) {
449 restartTime();
450 SimBlocParam(tt);
451 return 1;
452 }
453
454 if (iLast < 0) {
455 iLast = 0;
456 SimBlocReglage(tt);
457 return 1;
458 }
459
460 //if (stateChanged) { // emission d'un nouveau bloc reglage
461 // SimBlocReglage(tt);
462 // stateChanged=0;
463 // return 1;
464 //}
465
466
467 // Avons-nous un bloc bolo a transmettre ?
468 // duree contenue dans un bloc * numero de bloc...
469 t = numblock * per_echant() * (double)nb_per_block*2.;
470 if (iLast == 0 && ((TickCount() - tickStart)/60. < t/overSpeed)) return 0;
471
472 if ((TickCount() - tickStart)/60. > t/overSpeed + 10) {
473 numblock = (TickCount() - tickStart)/60.*overSpeed
474 / (per_echant() * (double)nb_per_block*2.);
475 iLast = 0;
476 }
477
478 if (iLast == 0) {
479 numblock++;
480 // On construit un scan...
481 PerformScan();
482 }
483
484 // On decide des blocs a envoyer...
485
486//#define need(n) \
487// (((simReglage.vitesse[n])>0) && ((numblock % (simReglage.vitesse[n])) ==0))
488//#define need(n) \
489// (((vitesse[n])>0) && ((numblock % (vitesse[n])) ==0))
490#define need(n) \
491 ((emission_blocks(n))!=0)
492
493 if (iLast < 1 && need(block_param)) {
494 iLast = 1;
495 SimBlocParam(tt);
496 return 1;
497 }
498
499 if (iLast < 2 && need(block_reglage)) {
500 iLast = 2;
501 SimBlocReglage(tt);
502 return 1;
503 }
504
505 if (iLast < 4 && need(block_gps)) {
506 iLast = 4;
507 SimBlocGPS(tt);
508 return 1;
509 }
510
511 if (iLast < 5 && need(block_une_periode)) {
512 iLast = 5;
513 SimBloc1Per(tt);
514 return 1;
515 }
516
517 if (iLast < 6 && need(block_bolo)) {
518 iLast = 6;
519 SimBlocBolo(tt);
520 return 1;
521 }
522
523 if (iLast < 7 && need(block_bolo_comprime)) {
524 iLast = 7;
525 SimBlocBoloComp(tt);
526 return 1;
527 }
528
529 if (iLast < 8 && need(block_sst)) {
530 iLast = 8;
531 SimBlocSST(tt);
532 return 1;
533 }
534
535 if (iLast < 9 && need(block_sst_comprime)) {
536 iLast = 9;
537 SimBlocSSTComp(tt);
538 return 1;
539 }
540
541 iLast = 0;
542 return 0;
543}
544
545#define verif_bolo {bolo=tc[2];if(bolo>nb_max_bolo) {printf("erreur: telecommande pour bolo=%d >nb_max_bolo \n",bolo);return;}}
546
547void simul_telecommandes_reduites(unsigned char* tc)
548{
549 int bolo,a;//,i;
550 if(tc[0]!=tc_reduite) return;
551
552 stateChanged = 1;
553 switch(tc[1])
554 {
555 case tc2_bolo_dacV:
556 verif_bolo;a=dac_V(simReglage.bolo[bolo]);
557 a=new_val_dac(a,tc[3]);
558 simReglage.bolo[bolo].mot1=(simReglage.bolo[bolo].mot1&0xfff000ff) | (a<<8);
559 printf("Simul: nouvelle valeur dacV=%d \n",a);
560 break;
561
562 case tc2_bolo_dacI:
563 verif_bolo;a=dac_I(simReglage.bolo[bolo]);
564 a=new_val_dac(a,tc[3]);
565 simReglage.bolo[bolo].mot1=(simReglage.bolo[bolo].mot1&0x000fffff) | (a<<20);
566 printf("Simul: nouvelle valeur dacI=%d \n",a);
567 break;
568
569 case tc2_bolo_gain:
570 verif_bolo;
571 simReglage.bolo[bolo].mot1=(simReglage.bolo[bolo].mot1&0xffffff00) | tc[3];
572
573 printf("Simul: nouvelle valeur gain=%x \n",tc[3]);
574 break;
575
576 case tc2_bolo_dacT:
577 verif_bolo;a=dac_T(simReglage.bolo[bolo]);
578 a=new_val_dac(a,tc[3]);
579 simReglage.bolo[bolo].mot2=(simReglage.bolo[bolo].mot2&0xfff000ff) | (a<<8);
580 printf("Simul: nouvelle valeur dacT=%d \n",a);
581 break;
582
583 case tc2_bolo_dacL:
584 verif_bolo;a=dac_L(simReglage.bolo[bolo]);
585 a=new_val_dac(a,tc[3]);
586 simReglage.bolo[bolo].mot2=(simReglage.bolo[bolo].mot2&0x000fffff) | (a<<20);
587 printf("Simul: nouvelle valeur dacL=%d \n",a);
588 break;
589
590 case tc2_bolo_voie:
591 verif_bolo;
592 simReglage.bolo[bolo].mot2=(simReglage.bolo[bolo].mot2&0xffffff00) | tc[3];
593
594 printf("Simul: nouvelle valeur voie=%d \n",tc[3]);
595 break;
596 case tc2_horloge:
597 switch(tc[2])
598 {
599 case tc3_commande :
600 printf(" commande generale numero =%d \n",tc[3]);
601 //commande_directe(tc[3]);
602 break;
603 case tc3_periode :
604 simReglage.horloge.periode=tc[3];
605 printf("Simul: nouvelle valeur periode=%d \n",tc[3]);
606 break;
607 case tc3_nb_mesures : simReglage.horloge.nb_mesures=tc[3];
608 printf("Simul: nouvelle valeur nb_mesures=%d \n",tc[3]);
609 break;
610 case tc3_temp_mort : simReglage.horloge.temp_mort=tc[3];
611 printf("Simul: nouvelle valeur temp_mort=%d \n",tc[3]);
612 break;
613 case tc3_flag : simReglage.horloge.flag=tc[3];
614 printf("Simul: nouvelle valeur flag=%d \n",tc[3]);
615 break;
616 default :
617 // i=tc[2]-tc3_vitesse;
618 // if( (i>=0) && (i<nb_type_blocks) )
619 // simReglage.vitesse[i]=tc[3];
620 // printf("Simul: nouvelle valeur block%d vitesse %d \n",i,tc[3]);
621 break;
622 }
623 restartTime();
624 case tc2_regul:
625 {
626 char * pt;
627 printf("Simul: tc reduite: ecrit regul val=%d pour position %d \n",tc[3],tc[2]);
628 pt=(char*)(&(simReglage.regul[0]));
629 pt=pt+(tc[2]);
630 if(tc[2]<nombre_de_regul*sizeof(regul_bolo)) *pt=tc[3];
631 else printf("erreur telecommande reduite tc2_regul \n");
632 }
633 break;
634
635
636
637 case tc2_auto_bolo:
638 {
639 char * pt1;
640 char * pt2;
641 printf(" tc reduite: ecrit autobolo val=%d pour position %d \n",tc[3],tc[2]);
642 pt1=(char*)(&(simReglage.autom[0]));pt2=pt1;
643 pt1=pt1+(tc[2]); /* pointeur sur l'octet a ecrire */
644 pt2=pt2+(tc[2]&0x1c)+3; /* le pointeur sur le mode (dernier octet pour transputer) */
645 if(tc[2]<nombre_de_voies*sizeof(auto_bolo))
646 {
647 *pt1=tc[3];
648 *pt2=*pt2 | (char)(((long)tc[2]&0x3)<<4) | 0x80;
649 /* le code de commande dans les 4 bit de pd fort du mode */
650 /* plus un 0X80 pour dire qu'il y a une commande */
651 }
652 else printf("erreur telecommande reduite tc2_auto_bolo \n");
653
654 }
655 break;
656
657 case tc2_auto_dilu:
658 {
659 char * pt1;
660 char * pt2;
661 printf(" tc reduite: ecrit autodilu val=%d pour position %d \n",tc[3],tc[2]);
662 pt1=(char*)(&(simReglage.dilu));pt2=pt1;
663 pt1=pt1+(tc[2]); /* pointeur sur l'octet a ecrire */
664 if(tc[2]<sizeof(auto_dilu))
665 {
666 *pt1=tc[3] | 0x80;
667 /* plus un 0X80 pour dire qu'il y a une commande */
668 /* les valeurs des parametres d'auto_dilu ne depassent pas 128 */
669 }
670 else printf("erreur telecommande reduite tc2_auto_bolo \n");
671
672 if(pt1==(char*)(&(simReglage.dilu.transmission)) )
673 {
674 simReglage.dilu.transmission = simReglage.dilu.transmission & 0x7f;
675 }
676 }
677 break;
678
679 default : break;
680 }
681}
682
683// microvolts -> codage
684
685def_gains
686#define gainbrut(aa) ((char)(((aa).mot1&0x1f)))
687#define gain(aa) gains_reels[gainbrut(aa)]
688#define codsignal(muv) muv/1.e5 * (65536.*gain(simReglage.bolo[i]))
689
690void PerformScan()
691{
692 int nb_coups,aa;
693 int i,j;
694 float ra,dec,dra,ddec;
695 struct localgeom lastGeom;
696 struct localgeom nextGeom;
697 double perech;
698 double secondes;
699
700 nb_coups= simReglage.horloge.nb_mesures/2 - simReglage.horloge.temp_mort;
701 aa = (nb_coups<<14) + (nb_coups*190) ;
702 perech = per_echant();
703
704 for (j=0; j<nb_per_block*2; j++) { // mesures dans le bloc
705 secondes = (numblock * nb_per_block*2 + j) * perech;
706 for(i=0;i<5;i++) { // bolometre
707 double muv;
708 if (!fmisdat) {
709 muv = 5 * sin((numblock * nb_per_block*2 + j)/(100.*(1+i)));
710 } else {
711 double noise;
712 muv = lastData.bolo[i].power + (nextData.bolo[i].power-lastData.bolo[i].power)*iGensample/oversample;
713 muv *= (3.e8/bolfreq[i]) * (3.e8/bolfreq[i]);
714 noise = GauRnd(0, whiteNoise[i]/sqrt(2*perech));
715 muv *= 1e6; noise *= 1e6;
716 muv *= 3.e8 * 1.e3; noise *= 3.e8 * 1.e3;
717 // Amortissement avec valeur precedente, sans le bruit
718 muv += boloLast[i] * exp(-perech/boloTimeCst);
719 boloLast[i] = muv;
720 muv += noise;
721 }
722 // add noise
723 //muv += GauRnd(0, whiteNoise);
724
725 databolo[i][j] = ((int)(((codsignal(muv) * (1-2*((j+1)%2)))*nb_coups)) >> 1 )+aa;
726 }
727 calcLocalGeom(&lastData, &lastGeom);
728 calcLocalGeom(&nextData, &nextGeom);
729 dra = nextGeom.ra-lastGeom.ra;
730 ddec = nextGeom.dec-lastGeom.dec;
731
732 // renormalize dra, ddec to have a norm of 1 arc min on sky
733 lastGeom.dra = dra;
734 lastGeom.ddec = ddec;
735 if (lastGeom.dra>12) lastGeom.dra-=24;
736 if (lastGeom.dra<-12) lastGeom.dra+=24;
737 {
738 double nrm;
739 nrm = lastGeom.dra*15*cos(lastGeom.dec * 3.1415926/180);
740 nrm = sqrt(lastGeom.ddec*lastGeom.ddec + nrm*nrm)*60;
741 lastGeom.ddec /= nrm;
742 lastGeom.dra /= nrm;
743 }
744
745 ra = lastGeom.ra + dra*iGensample/oversample;
746 dec = lastGeom.dec + ddec*iGensample/oversample;
747 // if (iGensample == 25 && iBuf == 0 && !hasDebugInfo) {
748 // hasDebugInfo=ra*1000;
749 // }
750 /* if (fabs(ra-17.41980)<0.00002 && fabs(dec-85.95694)<0.00002) {
751 //strcpy(debugInfo, " Debug 1");
752 hasDebugInfo=iGensample;
753 }*/
754
755 SSTSignal(ra,dec,lastGeom.dra,lastGeom.ddec,lastGeom.ora,lastGeom.odec,datadiod[j],secondes);
756
757 latitude = lastData.pos.lat + (nextData.pos.lat-lastData.pos.lat)*iGensample/oversample;
758 longitude = lastData.pos.lon + (nextData.pos.lon-lastData.pos.lon)*iGensample/oversample;
759 mjd = lastData.mjd + (nextData.mjd-lastData.mjd)*iGensample/oversample;
760 //JD+YGH on change le tps entre 2 echantillons pour simuler une periode differente
761 //mjd = lastData.mjd + (nextData.mjd-lastData.mjd)*iGensample/oversample/1.5;
762
763
764 iGensample ++;
765 if (iGensample >= oversample) {
766 iGensample = 0;
767 iBuf++;
768 if (iBuf >= szmisdatbuf-1) {
769 needBuf = curBuf;
770 curBuf = 1-curBuf;
771 iBuf = 0;
772 }
773 }
774 }
775}
776
777
778void SimBlocParam(tmtc* tt)
779{
780
781 /*
782 simParam.n_max_bolo = nb_max_bolo;
783 simParam.n_per_block = nb_per_block;
784 simParam.n_max_mes_per = nb_max_mes_per;
785
786 simParam.nb_bolo = 6;
787
788 for (i=0; i<nb_max_bolo; i++)
789 simParam.bolo[i].bolo_code_util = (i<6) ? bolo_normal_transmis : bolo_hors_service;
790
791 strcpy(simParam.bolo[0].bolo_nom, "bolo 1");
792 strcpy(simParam.bolo[1].bolo_nom, "bolo 2");
793 strcpy(simParam.bolo[2].bolo_nom, "bolo 3");
794 strcpy(simParam.bolo[3].bolo_nom, "bolo 4");
795 strcpy(simParam.bolo[4].bolo_nom, "bolo 5");
796 strcpy(simParam.bolo[5].bolo_nom, "bolo 6");
797 for (i=6; i<nb_max_bolo; i++)
798 simParam.bolo[i].bolo_nom[0] = 0;
799
800
801 for (i=0; i<nb_max_bolo; i++)
802 simParam.bolo[i].bolo_bebo = 10;
803
804 */
805 ((block_type_param*)(&tt->vi.btt))->param = simParam;
806 valide_block(&tt->vi.btt,block_param,numblock);
807}
808
809void SimBlocReglage(tmtc* tt)
810{
811 ((block_type_reglage*)(&tt->vi.btt))->reglage = simReglage;
812 valide_block(&tt->vi.btt,block_reglage,numblock);
813}
814
815
816
817
818
819
820void SimBlocBolo(tmtc* tt)
821{
822 int i,j;
823 block_type_bolo* blk = (block_type_bolo*)(&tt->vi.btt);
824
825 for (j=0; j<nb_per_block*2; j++) { // mesures dans le bloc
826 for(i=0;i<5;i++) { // bolometre
827 blk->data_bolo[i][j] = databolo[i][j];
828 }
829 }
830
831 for(i=5;i<nb_bolo_util;i++)
832 for (j=0; j<nb_per_block*2; j++) {
833 blk->data_bolo[i][j] = 0;
834 }
835
836 valide_block(&tt->vi.btt,block_bolo,numblock);
837}
838
839void SimBlocBoloComp(tmtc* tt)
840{
841 block_type_bolo_comprime blk2;
842 block_type_bolo* blk;
843 int i;
844 SimBlocBolo(tt);
845 blk = (block_type_bolo*) &(tt->vi.btt);
846 for(i=0;i<nb_bolo_util;i++) {
847 compress_7_2((unsigned long*)blk->data_bolo[i],(unsigned long*)blk2.data_bolo[i],nb_per_block*2,1);
848 }
849 valide_block((block_type_modele*)&blk2,block_bolo_comprime,numero_block(blk));
850 memcpy(&tt->vi.btt,&blk2,sizeof(blk2));
851}
852
853void SimBloc1Per(tmtc* tt)
854{
855 int i,j;
856
857 block_type_une_periode* blk = (block_type_une_periode*)(&tt->vi.btt);
858
859 for (i=0; i<6; i++)
860 for (j=0; j<simReglage.horloge.nb_mesures; j++) {
861 short value = (j*2<simReglage.horloge.nb_mesures) ?
862 4000 + 2000*exp(-j) : -4000-2000*exp(-(j-simReglage.horloge.nb_mesures/2));
863 blk->bol_per[i][j] = (value + 0x4000) ^ 0x7fff;
864 }
865 for (i=6; i<nb_max_bolo; i++)
866 for (j=0; j<simReglage.horloge.nb_mesures; j++)
867 blk->bol_per[i][j] = 0;
868 valide_block(&tt->vi.btt,block_une_periode,numblock);
869}
870
871// Format bloc GPS
872// $GPGGA,hhmmss.ss,ddmm.mmmm,n,dddmm.mmmm,e,q,ss,y.y,a.a,z,
873
874void SimBlocGPS(tmtc* tt)
875{
876 int h,m,s;
877 int dlo, dla;
878 double mlo, mla, j;
879 char clo, cla;
880 block_type_gps* blk = (block_type_gps*)(&tt->vi.btt);
881 j = (mjd - (int)mjd + 0.5) * 24;
882 h = j;
883 j = (j-h)*60;
884 m = j;
885 s = (j-m)*60;
886 cla = 'N'; j = latitude;
887 if (latitude < 0) {
888 cla = 'S'; j=-latitude;
889 }
890 dla = j;
891 mla = (j-dla)*60;
892 clo = 'E'; j = longitude;
893 if (longitude < 0) {
894 clo = 'W'; j=-longitude;
895 }
896 dlo = j;
897 mlo = (j-dlo)*60;
898 //strcpy(blk->gps, "$GPGGA,042232,3827.7653,N,00134.2222,E,1,07,2.3,32310.3,M\n");
899 //sprintf(blk->gps, "$GPGGA,%02d%02d%02d,%02d%07.4f,%1c,%02d%07.4f,%1c,1,07,2.3,32310.3,M\n",
900 sprintf(blk->gps, "$%02d%02d%02d,%02d%07.4f,%1c,%02d%07.4f,%1c,1,07,02.3,000,32310.3,M,32280,M,,\n",
901 h,m,s, dla, mla, cla, dlo, mlo, clo);
902
903 valide_block(&tt->vi.btt,block_gps,numblock);
904}
905
906void code_sst(block_type_sst* blk, int i, int* diodes); // diodes = tableau a 48 entrees
907
908// diodpermut[i] = channel de la diode i
909static int diodpermut[46]=
910 { 8,24,40, 9,25,41,10,26,42,11,
911 27,43,16,32, 1,17,33, 2,18,34,
912 3,19,35,12,28,44,13,29,45,14,
913 30,46,15,31,47,20,36, 5,21,37,
914 6,22,38, 7,23,39};
915 // voies 0 et 4 non connectees, voie 1 en panne.
916
917void code_sst(block_type_sst* blk, int i, int* diodes) {
918 int j; // 0-5 : numero du bloc de 8 diodes
919 int k; // 0-2 : indice du bloc de 4 bits (une diode = 12 bits = 3 blocs de 4 bits), MSB=0
920 int l; // 0-7 : indice de la diode dans son bloc (8 diodes * 4 bits = 1 mot de 32 bits)
921 int dd[46];
922
923 // numero de la diode (0-47) = j*8+l;
924 // indice dans le bloc sst du mot de 32 bits (0-17) = j*3+k;
925 // indice dans mot de 32 bits du premier bit utile = 4*l;
926
927 // Permutation des diodes
928 for (j=0; j<46; j++) {
929 dd[j] = diodes[j];
930 diodes[j]=0;
931 }
932 diodes[46] = diodes[47] = 0;
933
934 for (j=0; j<46; j++) {
935 diodes[diodpermut[j]] = dd[j];
936 }
937 diodes[0] = 23 + 2048; // canal non connecte
938 diodes[4] = 45 + 2048; // canal non connecte
939 diodes[1] = 1234 + 2048; // canal 1 en panne pour le moment.
940
941 for (j=0; j<6; j++)
942 for (k=0; k<3; k++) {
943 blk->sst[i][j*3+k] = 0;
944 for (l=0; l<8; l++) {
945 short bit4 = (diodes[j*8+l] >> 4*(2-k)) & 0xF;
946 blk->sst[i][j*3+k] += (bit4 << (4*l));
947 }
948 }
949}
950
951
952void SimBlocSST(tmtc* tt)
953{
954 int i;
955 block_type_sst* blk = (block_type_sst*)(&tt->vi.btt);
956
957 for (i=0; i<nb_per_block*2; i++) {
958 code_sst(blk, i, datadiod[i]);
959 }
960
961 valide_block(&tt->vi.btt,block_sst,numblock);
962}
963
964#define place_paquet(i,j) ((i/8) * 24 + j*8 + (i%8) )
965
966void SimBlocSSTComp(tmtc* tt)
967{
968 block_type_sst_comprime blk2;
969 block_type_sst* blk;
970 int j,k,jc;
971 unsigned long sst_vrai[nb_per_block*2];
972 unsigned long a,b0,b1,b2;
973 SimBlocSST(tt);
974 blk = (block_type_sst*) &(tt->vi.btt);
975/* pour nb_per_block=36 periodes completes on a 72 points a comprimer pour chaque diode */
976
977 jc=0;
978 for(j=0;j<48;j++) /* jc = bolo_comprime j=bolo normal */
979 {
980 if( (j!=0) && (j!=4) )
981 {
982 for(k=0;k<nb_per_block*2;k++) /* boucle sur les demi entières */
983 {
984 a=place_paquet(j,0);
985 b0= ( blk->sst[k][a/8] >>( (a%8)*4) ) & 0xf;
986 a=place_paquet(j,1);
987 b1= ( blk->sst[k][a/8] >>( (a%8)*4) ) & 0xf;
988 a=place_paquet(j,2);
989 b2= ( blk->sst[k][a/8] >>( (a%8)*4) ) & 0xf;
990
991 sst_vrai[k]=( (b0<<8) | (b1<<4) | b2 ) ;
992 }
993 compress_4_1(sst_vrai,blk2.sst[jc],nb_per_block*2,1);
994 jc++;
995 }
996 }
997
998 valide_block((block_type_modele*)&blk2,block_sst_comprime,numero_block(blk));
999 memcpy(&tt->vi.btt,&blk2,sizeof(blk2));
1000}
1001
1002
1003#define SIDRATE 0.9972695677
1004#define refTS (13*3600. + 10.*60. + 33.1)
1005#include "aa_hadec.h"
1006
1007void calcLocalGeom(struct misdat* md, struct localgeom* geom)
1008{
1009 // Calcul temps sideral local
1010 // GMT TS Pour JD = 2450000 = 13h10m33.1s
1011 double ts; // secondes
1012 double ha, dec, lat, alt, az;
1013
1014 // La longitude est vers l'EST
1015
1016 ts = md->mjd * 86400. / SIDRATE + refTS + (md->pos.lon/15. * 3600.);
1017 ts = ts - 86400.*(int)(ts/86400);
1018
1019 // Elevation = 41 degres... -> Ra & Dec
1020 alt = 41. * M_PI/180;
1021 lat = md->pos.lat * M_PI/180;
1022 az = md->pos.azimut * M_PI/180;
1023
1024 aa_hadec(lat, alt, az, &ha, &dec);
1025
1026 geom->ra = - (ha * 180. / M_PI / 15) + (ts/3600.);
1027 geom->dec = (dec * 180. / M_PI);
1028
1029 alt = (41. + 1./60.) * M_PI/180;
1030
1031 aa_hadec(lat, alt, az, &ha, &dec);
1032
1033 geom->ora = - (ha * 180. / M_PI / 15) + (ts/3600.) - geom->ra;
1034 geom->odec = (dec * 180. / M_PI) - geom->dec;
1035// en fait, prendre plus tard l'orthogonal du vecteur de balayage, dans le
1036// sens le plus proche de cette valeur (indetermination sens de rotation)
1037// pour tenir compte de la pendulation
1038
1039}
1040
1041
1042// Etalonnage approx : 3e8 V/W
1043// Cte temps 10 ms
1044
1045#define DeuxPi (2* M_PI)
1046
1047
1048// Generation aleatoire gaussienne de centre "am" et de sigma "s".
1049
1050double GauRnd(double am, double s)
1051{
1052double x,A,B;
1053
1054LAB10:
1055A = drand01();
1056if ( A == 0. ) goto LAB10;
1057B = drand01();
1058x = am + s * sqrt(-2.*log(A))*cos(DeuxPi*B);
1059return(x);
1060}
Note: See TracBrowser for help on using the repository browser.