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

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

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