source: Sophya/trunk/Poubelle/archediab.old/archediab.sources/c/recons_sst.c@ 649

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

archediab 28 vol

File size: 15.3 KB
RevLine 
[649]1#include "diabolo.h"
2#include "senseur_stellaire.h"
3#include "recons_sst.h"
4#include "FindPeri.h"
5#include "stabilite.h"
6FILE* f;
7long nstarcatal,nposaxesimu;
8struct star {
9 float ra; // en heures dans le catalogue de depart
10 float dec;// en degres dans le catalogue de depart
11 short mag; // mag * 100
12};
13#define nstpturn 100
14struct sstgraph {
15 int netdet; // nombre d'etoiles detectees
16 float raax; // coordonnees equatoriales de l'axe de rotation du telescope
17 float decax;
18 long numcat[nstpturn]; //numero de l'etoile dans le fichier de depart
19 float phase[nstpturn]; //phase de l'etoile (angle/Nord - Est=+90¡ - ds plan perp axe rot ballon)
20 float cosinus[nstpturn]; // de l'angle entre la direction de l'etoile et l'axe de rotation du ballon
21 float fluxcat[nstpturn]; //flux de l'etoile en unites arbitraires
22 int icl[nstpturn]; //ordre de l'etoile pour un classement par phases croissantes
23};
24struct star* stars;
25struct sstgraph* listord;
26int fichierslus=0;
27float instant[nmemechant],sfluxvu[nmemechant],barette[nmemechant],phd0[nmemechant];
28float deltap[nmemechant],rapflx[nmemechant];
29float deltapcat[nstpturn],rapflxcat[nstpturn];
30float cases[ncasdiftps];float maxcas;float somcas;
31float sflxv,maintenant,laplushaute;
32float tdeb;int flagstar;int flagwind;
33float separtps,separenz;
34float Period=10.,NewPeriod;float Phase;int Sens=9999;
35float WidthCas,StepCas;
36float tpsminfen,tpsmaxfen;
37float Liminffen,Limsupfen;
38double CoupFlux;
39float Papa,poidsvois;
40float tauxper;long nstastat,ndsper;
41int flagChPer,flagvuefiltree;
42int nptsfenper,nptsfenphs;
43float somtps,somecartps,somecarz;
44float memtemps,membar;
45int ncoups;int statcalcphase[ncasphase];int CasesVides;
46float pasenphase,ecarmoytps[ncasphase],ecarmoyz[ncasphase];
47float sini[ncasphase],cosini[ncasphase];
48float phaseret=-9999;
49float pluvieux;int nlistevue,numaccro,flcatmin;
50#define D2R .0174533 //1 degre = D2R radians
51#define H2R .261799 // 1 heure siderale = H2R radians
52#define P2 1.570796 // pi sur 2
53#define PI 3.14159265359 // pi
54#define DPI 6.2831853 // 2*pi
55//variables geometriques du senseur stellaire :
56
57#define sommetdeg 49. // 1/2 angle d'ouverture du cone moye balayŽ par la barette (en degres)
58#define ouvertdeg 1.44 // extension de la barette en delta (en degres)
59float sommet,ouvert; // memes angles en radians
60 float decaxe;
61 float raaxe;
62
63
64void init_recons_sst()
65
66//appele au debut par diabolo.c, puis a chaque reset de la fenetre findperiod
67{
68int i;
69double valeur=0;
70//printf("periode echantillonnage %f",gg->periode_echantillonage);
71pasenphase=360./(float)(ncasphase);
72for(i=0;i<nmemechant;i++)
73{
74instant[i]=0.;
75sfluxvu[i]=0.;
76barette[i]=0.;
77}
78for(i=0;i<ncasdiftps;i++)
79{
80cases[i]=0;
81}
82flagvuefiltree=0;
83//exec_FindPeriod(fenetre_Find_Period,ouverture,valeur);
84flagwind=0;
85maxcas=0;
86Period=LowSeek;Phase=-9999.;Sens=-9999;
87somcas=0;
88tdeb=0;
89flagstar=0;
90CoupFlux=10.;
91if(fenetre(fenetre_Find_Period)) CoupFlux=Rd_CpFl(fenetre_Find_Period);
92//printf("coupure en flux %d \n",(int)(CoupFlux));
93tauxper=0.;
94nstastat=0;
95ndsper=0;
96nptsfenper=0;
97nptsfenphs=0;
98memtemps=-1.;
99ncoups=0;
100somtps=0;somecartps=0;somecarz=0;
101
102 Liminffen=LowSeek;
103 Limsupfen=HighSeek;
104 Papa=10.;
105 if(fenetre(fenetre_Find_Period))
106 {
107 if(DonnePeriod() != -1)
108 {
109 Liminffen=DonnePeriod()-1;
110 Limsupfen=DonnePeriod()+1;
111 Papa=1.;
112 }
113 flagwind=1;
114 }
115
116 passargu(Liminffen,Limsupfen,Papa);
117 //nouveauD(fenetre_Find_Period,FindPeriod_id,"FindPeriod",exec_FindPeriod);
118 //nouveauD(fenetre_Phase_SST,phase_id,"phase",exec_phase);
119 //corFC flagwind=0;
120 raz_fen_phase();
121 if (fenetre(fenetre_Find_Period)) {
122 selectgra(fenetre_Find_Period);
123 graph->xmin=Liminffen;
124 graph->xmax=Limsupfen;
125 graph->xpas=Papa;
126 retrace(fenetre_Find_Period);
127 ecritD(fenetre_Find_Period,FP_vperesti," %7.3f",DonnePeriod());
128 }
129
130 tpsminfen=Liminffen;
131 tpsmaxfen=Limsupfen;
132StepCas=(HighSeek-LowSeek)/(float)(ncasdiftps);
133WidthCas=2*StepCas;
134tabulesincos();
135CasesVides=ncasphase;
136//lecture des catalogues STELLAIRES!
137if(fichierslus==0)
138{
139fichierslus=1;
140printf(" \n");
141 // Catalogue Dominique Yvon
142 {
143 short ref;
144 long n;
145 if (FSOpen("\phyp6HADY.dat",0,&ref)) {
146 printf("***ERREUR*** ouverture hyp6Hammam.dat\n");
147 return;
148 }
149 n = sizeof(long);
150 FSRead(ref, &n, (char*) &nstarcatal);
151 printf("nstarcatal = %d\n", nstarcatal);
152 n = nstarcatal*sizeof(struct star);
153 stars = (struct star*) malloc(n);
154 FSRead(ref, &n, (char*) stars);
155 FSClose(ref);
156 printf("Hyp6DY m<6 : read %d stars\n", nstarcatal);
157 }
158 for(i=0;i<5;i++)
159 {
160 printf("alpha %f,delta %f, mag %d \n",stars[i].ra,stars[i].dec,stars[i].mag);
161 }
162 {
163 short ref;
164 long n;
165 if (FSOpen("\psstgraphs.dat",0,&ref)) {
166 printf("***ERREUR*** ouverture sstgraphs.dat\n");
167 return;
168 }
169 n = sizeof(long);
170 FSRead(ref, &n, (char*) &nposaxesimu);
171 printf("nstarcatal = %d\n", nposaxesimu);
172 n = nposaxesimu*sizeof(struct sstgraph);
173 listord = (struct sstgraph*) malloc(n);
174 FSRead(ref, &n, (char*) listord);
175 FSClose(ref);
176 printf("sstgraph.dat : read %d positions simulees\n", nposaxesimu);
177 }
178 printf("nombre d'etoiles selectionnees %d\n",listord[0].netdet);
179 numaccro=0;flcatmin=1;
180}
181}
182
183void exec_recons_sst()
184{
185int i,icour;
186int cascour,casper,vudsper;
187int valid;
188float delbar,delebarre;
189float deltan,deltintin;
190float delcas,pocac;
191sflxv=0.;
192if(nfoundstars>0)
193{
194if(nfoundstars>MAXFOUNDSTARS) nfoundstars=MAXFOUNDSTARS;
195laplushaute=0.;
196for(i=0;i<nfoundstars;i++)
197{
198if(mfoundstars[i]>sflxv)
199{sflxv=mfoundstars[i];
200laplushaute=zfoundstars[i];
201}
202//printf("%d %f %f %f\n",i,tfoundstars[i],zfoundstars[i],mfoundstars[i]);
203}
204
205if(sflxv>CoupFlux)
206{
207maintenant=tfoundstars[0];
208membar=laplushaute;
209if(flagstar==0)
210 {
211 tdeb=maintenant;
212 flagstar=1;
213 }
214 if(maintenant-tdeb>LowSeek)
215 {
216 vudsper=0;
217 valid=0;
218 delbar=0;
219 delebarre=10000;
220 deltan=0.;
221 deltintin=1.;
222for(i=0;i<nmemechant;i++)
223{
224if(sfluxvu[i]==0.) break;
225separtps=maintenant-instant[i];
226separenz=laplushaute-barette[i];
227if(separtps>LowSeek && fabs(sfluxvu[i]-sflxv)/sflxv<.5)
228 {
229 if(separtps>tpsminfen && separtps<tpsmaxfen)
230 {
231 trace_point(fenetre_Find_Period);
232 nptsfenper++;
233 if(nptsfenper==QuorumRafrai && fenetre(fenetre_Find_Period))
234 { selectgra(fenetre_Find_Period);efface(fenetre_Find_Period);nptsfenper=0;}
235 }
236 // recherche de periode
237 if(separtps<HighSeek)
238 {
239 casper=(int)((Period-LowSeek)/StepCas);
240 cascour=(int)((separtps-LowSeek)/StepCas);
241 delcas=separtps-StepCas*(cascour+.5)-LowSeek;
242 pocac=fabs(delcas/StepCas);
243// printf("separtps %f delcas %f pocac %f\n",separtps,delcas,pocac);
244
245
246 if(cascour>0 && cascour<ncasdiftps-1)
247 {
248 cases[cascour]+=25./separtps*(1.-pocac);
249 if(delcas>0.){cases[cascour+1]+=25./separtps*pocac;}
250 else {cases[cascour-1]+=25./separtps*pocac;}
251 somcas+=25./separtps;
252 if(fabs(cascour-casper)<2)
253 {vudsper++;
254 }
255 if(cases[cascour]>maxcas)
256 {
257 maxcas=cases[cascour];
258 poidsvois=cases[cascour-1]+cases[cascour]+cases[cascour+1];
259 NewPeriod=LowSeek+
260 (cascour+.5+(cases[cascour+1]-cases[cascour-1])/poidsvois)*StepCas;
261 if(poidsvois>5*somcas/ncasdiftps){
262 flagChPer=0;
263 if(fabs(Period-NewPeriod)>1.)
264 {nstastat=0;ndsper=0;tauxper=0.;
265 { if(DonnePeriod() != -1 && fenetre(fenetre_Find_Period))
266 {selectgra(fenetre_Find_Period);efface(fenetre_Find_Period);nptsfenper=0;}}}
267 if(fabs(Period-NewPeriod)>.1) flagChPer=1;
268 if(fabs(Period-NewPeriod)>.010)
269 {raz_fen_phase();}
270 Period=NewPeriod;
271 }
272// printf("maxcas %f somcas/ncasdiftps %f Periode : %f\n",maxcas,
273// somcas/(float)(ncasdiftps),Period);
274 }
275 }
276 }
277 delbar=separenz;
278 deltan=separtps-Period;
279 if(DonnePeriod() != -1
280 && fabs(deltan-Period)<FourchTps && fabs(deltan-Period-2*deltintin)<FourchTps
281 && fabs(delebarre-2*delbar)<8 && valid==1)
282 {
283 tracephase();
284// membar=.3*barette[i]+.7*membar;
285 nptsfenphs++;
286 if(nptsfenphs==RafraiPhase) {raz_fen_phase();}
287// printf("periode %8.3f phases tpsdiodes %8.1f sens %d\n",DonnePeriod(),DonnePhase(),DonneSens());
288
289 }
290 if(DonnePeriod() != -1 && fabs(deltan)<FourchTps && fabs(delbar)<8)
291 {valid=1;delebarre=delbar;deltintin=deltan;}
292
293 }
294}
295nstastat++;
296if(vudsper!=0)ndsper++;
297if(nstastat>StatValidePeriode)
298{
299 tauxper=(float)(ndsper)/(float)(nstastat);
300
301// printf("efficacite periode %f\n",tauxper);
302}
303
304if(nstastat==50) {nstastat=0;ndsper=0;}
305if(DonnePeriod()==-1){
306raz_fen_phase();
307if(fenetre(fenetre_Find_Period) && tpsminfen != 10.){
308 selectgra(fenetre_Find_Period);
309 graph->xmin=10.;
310 graph->xmax=40.;
311 graph->xpas=10.;
312 Liminffen=graph->xmin;
313 Limsupfen=graph->xmax;
314 tpsminfen=graph->xmin;
315 tpsmaxfen=graph->xmax;
316 passargu(Liminffen,Limsupfen,Papa);
317 efface(fenetre_Find_Period);
318 retrace(fenetre_Find_Period);
319 nptsfenper=0;
320}
321
322} else
323{
324if(flagChPer==1){
325 selectgra(fenetre_Find_Period);
326 graph->xmin=DonnePeriod()-1.;
327 graph->xmax=DonnePeriod()+1.;
328 graph->xpas=1.;
329 Liminffen=graph->xmin;
330 Limsupfen=graph->xmax;
331 tpsminfen=graph->xmin;
332 tpsmaxfen=graph->xmax;
333 passargu(Liminffen,Limsupfen,Papa);
334 efface(fenetre_Find_Period);
335 retrace(fenetre_Find_Period);
336 nptsfenper=0;
337}
338}
339if(fenetre(fenetre_Find_Period))
340{selectgra(fenetre_Find_Period);ecritD(fenetre_Find_Period,FP_vperesti," %7.3f",DonnePeriod());}
341if(Liminffen==LowSeek && DonnePeriod() != -1 && fenetre(fenetre_Find_Period))
342{
343 selectgra(fenetre_Find_Period);
344 graph->xmin=DonnePeriod()-1.;
345 graph->xmax=DonnePeriod()+1.;
346 graph->xpas=1.;
347 Liminffen=graph->xmin;
348 Limsupfen=graph->xmax;
349 tpsminfen=graph->xmin;
350 tpsmaxfen=graph->xmax;
351 passargu(Liminffen,Limsupfen,Papa);
352 efface(fenetre_Find_Period);
353 retrace(fenetre_Find_Period);
354
355}
356}
357
358icour=0;
359pluvieux=maintenant;
360for(i=nmemechant-1;i>0;i--)
361 {
362 icour++;
363 if(sfluxvu[i-1]==0.) continue;
364 if(instant[i]<pluvieux)pluvieux=instant[i];
365 sfluxvu[i]=sfluxvu[i-1];
366 instant[i]=instant[i-1];
367 barette[i]=barette[i-1];
368 phd0[i]=phd0[i-1];
369 if(DonnePeriod()!=-1){if(maintenant-instant[i]>3.*DonnePeriod())
370 {
371 sfluxvu[i]=0.;
372 instant[i]=0.;
373 barette[i]=0.;
374 phd0[i]=0.;
375 }
376 }
377 }
378sfluxvu[0]=sflxv;
379instant[0]=maintenant;
380barette[0]=laplushaute;
381phd0[i]=0.;
382{float peri=DonnePeriod();
383if (peri>0.)phd0[i]=360./peri*(maintenant-peri*(int)(maintenant/peri));
384}
385// bloc recherche de l'axe ->identification des etoiles
386{float difftps=0.;
387
388if(DonnePeriod()>0. && (maintenant - pluvieux)>DonnePeriod())
389{
390nlistevue=0;
391for(i=0;i<nmemechant-1;i++)
392{
393 if(sfluxvu[i]==0. |maintenant-instant[i]>Period) break;
394difftps=(instant[i]-instant[i+1])*DPI/Period;
395nlistevue++;
396deltap[i]=difftps;
397rapflx[i]=sfluxvu[i+1]/sfluxvu[i];
398}
399//printf("%d etoiles detectees par tour\n",nlistevue);
400if(fenetre(fenetre_Find_Period))
401{selectgra(fenetre_Find_Period);ecritD(fenetre_Find_Period,Find_net," %d",nlistevue);}
402numaccro=1;
403if(numaccro==0) chercheaccro(flcatmin);
404
405
406
407
408}}
409// fin du bloc recherche de l'axe ->identification des etoiles
410
411}
412}
413}
414void chercheaccro(int seuil)
415{
416/*
417#define nstpturn 100
418struct sstgraph {
419 int netdet; // nombre d'etoiles detectees
420 float raax; // coordonnees equatoriales de l'axe de rotation du telescope
421 float decax;
422 long numcat[nstpturn]; //numero de l'etoile dans le fichier de depart
423 float phase[nstpturn]; //phase de l'etoile (angle/Nord - Est=+90¡ - ds plan perp axe rot ballon)
424 float cosinus[nstpturn]; // de l'angle entre la direction de l'etoile et l'axe de rotation du ballon
425 float fluxcat[nstpturn]; //flux de l'etoile en unites arbitraires
426 int icl[nstpturn]; //ordre de l'etoile pour un classement par phases croissantes
427};
428*/
429int i,j;int nstartour;float temps[nstpturn],fluxion[nstpturn];
430nstartour=0;
431for(i=0;i<nposaxesimu;i++)
432 {
433 for(j=0;j<listord[i].netdet-1;j++)
434 {
435 deltapcat[j]=listord[i].phase[j+1]-listord[i].phase[j];
436 rapflxcat[j]=listord[i].fluxcat[j+1]/listord[i].fluxcat[j];
437 }
438 }
439
440
441}
442void trace_point(int fen)
443{
444if(fenetre(fen)) symbole(fen, separtps, sflxv, 2, rond, 0, bleu);
445}
446void tracephase(void)
447{
448float tempsplie,tempsmoy;float peri;float ecartps,ecarz;
449peri=DonnePeriod();
450if (peri>0.)
451{
452tempsplie=360./peri*(maintenant-peri*(int)(maintenant/peri));
453if(!flagvuefiltree)
454{
455if(fenetre(fenetre_Phase_SST))
456{symbole(fenetre_Phase_SST,tempsplie,separtps-2.*peri, 10, carre, 0, bleu);
457symbole(fenetre_Phase_SST,tempsplie,.005*(separenz), 10, croix, 0, rouge);}
458}
459
460{
461if(ncoups==0)memtemps=((int)(tempsplie)/pasenphase)*pasenphase;
462 if(tempsplie-memtemps<pasenphase && tempsplie>=memtemps)
463 {
464 ncoups++;
465 somtps+=tempsplie;
466 somecartps+=separtps-2.*peri;
467 somecarz+=(separenz)*.005;
468 } else
469 {
470 if(ncoups>0)
471 {
472 tempsmoy=somtps/(float)(ncoups);
473// printf("temps,ecart : %f %f\n",tempsmoy,somecartps/ncoups);
474 ecartps =somecartps/(float)(ncoups);
475 ecarz = somecarz/(float)(ncoups);
476 charge_calc_phase(tempsmoy,ecartps,ecarz);
477if(fenetre(fenetre_Phase_SST) && flagvuefiltree)
478{ symbole(fenetre_Phase_SST,tempsmoy,ecartps, 10, carre, 0, bleu);
479 symbole(fenetre_Phase_SST,tempsmoy,ecarz, 10, croix, 0, rouge);
480}
481 ncoups=0;
482 somtps=0;
483 somecartps=0;
484 somecarz=0;
485 memtemps=((int)(tempsplie)/pasenphase)*pasenphase;
486 {
487 ncoups++;
488 somtps+=tempsplie;
489 somecartps+=separtps-2.*peri;
490 somecarz+=(separenz)*.005;
491 }
492 }
493 }
494
495}
496}}
497float DonnePeriod(void)
498{
499if(tauxper>.4){
500return Period;}
501else return -1.;
502}
503int CpFl(void)
504{
505return (int)(CoupFlux);
506}
507void raz_fen_phase()
508{
509int i;
510if (fenetre(fenetre_Phase_SST)){selectgra(fenetre_Phase_SST);
511efface(fenetre_Phase_SST);}
512memtemps=-1.;ncoups=0;nptsfenphs=0;phaseret=-9999.;Sens=-9999;
513CasesVides=ncasphase;for(i=0;i<ncasphase;i++) statcalcphase[i]=0;
514
515}
516void charge_calc_phase(float tempsmoy,float ecartps,float ecarz)
517{
518int caserne;
519caserne=(tempsmoy+pasenphase*.5)/pasenphase;
520if(caserne==ncasphase)caserne=0;
521if(caserne < ncasphase && caserne >= 0)
522{
523if(statcalcphase[caserne]<5 | (fabs(ecarmoytps[caserne]-ecartps)<.02 &&
524fabs(ecarmoyz[caserne]-ecarz)<3.))
525{ecarmoytps[caserne]*=statcalcphase[caserne];
526ecarmoyz[caserne]*=statcalcphase[caserne];
527ecarmoytps[caserne]+=ecartps;
528ecarmoyz[caserne]+=ecarz;
529statcalcphase[caserne]++;
530ecarmoytps[caserne]/=statcalcphase[caserne];
531ecarmoyz[caserne]/=statcalcphase[caserne];
532CasesVides=calc_phase();
533}}}
534void tabulesincos(void)
535{
536int i;float pasenrad;
537pasenrad=pasenphase*.01745329;
538for(i=0;i<ncasphase;i++)
539{
540sini[i]=sin(pasenrad*(i));
541cosini[i]=cos(pasenrad*(i));
542}
543}
544int calc_phase(void)
545{
546int i,vides;float cocoz,sisiz,cocotps,sisitps;
547float phasez,phasetps;float sensrot;
548static float phasemem;
549vides =0;
550for(i=0;i<ncasphase;i++)
551{
552if(statcalcphase[i]==0) vides++;
553}
554
555cocoz=0;
556sisiz=0;
557cocotps=0;
558sisitps=0;
559for(i=0;i<ncasphase;i++)
560{
561if(statcalcphase[i]>0)
562{
563cocoz+=cosini[i]*ecarmoyz[i];
564sisiz+=sini[i]*ecarmoyz[i];
565cocotps+=cosini[i]*ecarmoytps[i];
566sisitps+=sini[i]*ecarmoytps[i];
567}
568}
569phasez=-atan2(cocoz,sisiz);
570phasetps=-atan2(cocotps,sisitps);
571sensrot=phasez-phasetps;
572if(sensrot>3.14)sensrot-=6.28;
573if(sensrot<-3.14)sensrot+=6.28;
574if(sensrot>0.) Sens =1;
575if(sensrot<=0.) Sens =-1;
576//arbitraire, ajuste sur la simulation
577Phase=phasez*57.29577;
578phaseret=-9999.;
579if(CasesVides<3 && fabs(Phase-phasemem)<10.){ phaseret=Phase;}
580if(CasesVides<4) phasemem=Phase;
581//printf("phases : numdiodes %f, tpsdiodes %f sens %d\n",phasez,phasetps,Sens);
582return vides;}
583float DonnePhase(void)
584{
585return phaseret;
586}
587int DonneSens(void)
588{
589if(CasesVides<3){return Sens;}else {return -9999;}
590}
Note: See TracBrowser for help on using the repository browser.