#include "diabolo.h" 
#include "senseur_stellaire.h" 
#include "recons_sst.h" 
#include "FindPeri.h" 
#include "stabilite.h"
FILE* f;
long nstarcatal,nposaxesimu;
struct star {
  float ra; // en heures dans le catalogue de depart
  float dec;// en degres dans le catalogue de depart
  short mag; // mag * 100
};
#define nstpturn 100
struct sstgraph {
  int netdet; // nombre d'etoiles detectees
  float raax; // coordonnees equatoriales de l'axe de rotation du telescope
  float decax;
  long numcat[nstpturn]; //numero de l'etoile dans le fichier de depart
  float phase[nstpturn]; //phase de l'etoile (angle/Nord - Est=+90ˇ - ds plan perp axe rot ballon)
  float cosinus[nstpturn]; // de l'angle entre la direction de l'etoile et l'axe de rotation du ballon
  float fluxcat[nstpturn]; //flux de l'etoile en unites arbitraires
  int	icl[nstpturn]; //ordre de l'etoile pour un classement par phases croissantes
};
struct star* stars;
struct sstgraph* listord;
int fichierslus=0;
float instant[nmemechant],sfluxvu[nmemechant],barette[nmemechant],phd0[nmemechant];
float deltap[nmemechant],rapflx[nmemechant];
float deltapcat[nstpturn],rapflxcat[nstpturn];
float cases[ncasdiftps];float maxcas;float somcas;
float sflxv,maintenant,laplushaute;
float tdeb;int flagstar;int flagwind;
float separtps,separenz;
float Period=10.,NewPeriod;float Phase;int Sens=9999;
float WidthCas,StepCas;
float tpsminfen,tpsmaxfen;
float Liminffen,Limsupfen;
double CoupFlux;
float Papa,poidsvois;
float tauxper;long nstastat,ndsper;
int flagChPer,flagvuefiltree;
int nptsfenper,nptsfenphs;
float somtps,somecartps,somecarz;
float memtemps,membar;
int ncoups;int statcalcphase[ncasphase];int CasesVides;
float pasenphase,ecarmoytps[ncasphase],ecarmoyz[ncasphase];
float sini[ncasphase],cosini[ncasphase];
float phaseret=-9999;
float pluvieux;int nlistevue,numaccro,flcatmin;
#define D2R .0174533 //1 degre = D2R radians
#define H2R .261799 // 1 heure siderale = H2R radians
#define P2 1.570796 // pi sur 2
#define PI 3.14159265359 // pi
#define DPI 6.2831853 // 2*pi
//variables geometriques du senseur stellaire :

#define sommetdeg 49. // 1/2 angle d'ouverture du cone moye balayŽ par la barette (en degres)
#define ouvertdeg 1.44 // extension de la barette en delta (en degres)
float sommet,ouvert; // memes angles en radians
   float decaxe;
   float raaxe;


void init_recons_sst()

//appele au debut par diabolo.c, puis a chaque reset de la fenetre findperiod
{
int i;
double valeur=0;
//printf("periode echantillonnage %f",gg->periode_echantillonage);
pasenphase=360./(float)(ncasphase);
for(i=0;i<nmemechant;i++)
{
instant[i]=0.;
sfluxvu[i]=0.;
barette[i]=0.;
}
for(i=0;i<ncasdiftps;i++)
{
cases[i]=0;
}
flagvuefiltree=0;
//exec_FindPeriod(fenetre_Find_Period,ouverture,valeur);
flagwind=0;
maxcas=0;
Period=LowSeek;Phase=-9999.;Sens=-9999;
somcas=0;
tdeb=0;
flagstar=0;
CoupFlux=10.;
if(fenetre(fenetre_Find_Period)) CoupFlux=Rd_CpFl(fenetre_Find_Period);
//printf("coupure en flux %d \n",(int)(CoupFlux));
tauxper=0.;
nstastat=0;
ndsper=0;
nptsfenper=0;
nptsfenphs=0;
memtemps=-1.;
ncoups=0;
somtps=0;somecartps=0;somecarz=0;

				Liminffen=LowSeek;
				Limsupfen=HighSeek;
				Papa=10.;
				if(fenetre(fenetre_Find_Period))
				 {
				 if(DonnePeriod() != -1)
				  {
				  Liminffen=DonnePeriod()-1;
				  Limsupfen=DonnePeriod()+1;
				  Papa=1.;
				  }
				 flagwind=1;
				 }
				
				passargu(Liminffen,Limsupfen,Papa);
					//nouveauD(fenetre_Find_Period,FindPeriod_id,"FindPeriod",exec_FindPeriod);
					//nouveauD(fenetre_Phase_SST,phase_id,"phase",exec_phase);
				//corFC	flagwind=0;
				raz_fen_phase();
 				if (fenetre(fenetre_Find_Period)) {
   					selectgra(fenetre_Find_Period);
    				graph->xmin=Liminffen;
    				graph->xmax=Limsupfen;
    				graph->xpas=Papa;
    				retrace(fenetre_Find_Period);
					ecritD(fenetre_Find_Period,FP_vperesti," %7.3f",DonnePeriod());
    			}

    tpsminfen=Liminffen;
    tpsmaxfen=Limsupfen; 
StepCas=(HighSeek-LowSeek)/(float)(ncasdiftps);
WidthCas=2*StepCas;
tabulesincos();
CasesVides=ncasphase;
//lecture des catalogues STELLAIRES!
if(fichierslus==0)
{
fichierslus=1;
printf(" \n");
  // Catalogue Dominique Yvon
  {
     short ref;
     long n;
     if (FSOpen("\phyp6HADY.dat",0,&ref)) {
         printf("***ERREUR*** ouverture hyp6Hammam.dat\n");
         return;
     }
     n = sizeof(long);
     FSRead(ref, &n, (char*) &nstarcatal);
   printf("nstarcatal = %d\n", nstarcatal);
     n = nstarcatal*sizeof(struct star);
     stars = (struct star*) malloc(n);
     FSRead(ref, &n, (char*) stars);
     FSClose(ref);
     printf("Hyp6DY m<6 : read %d stars\n", nstarcatal);
  }
   for(i=0;i<5;i++)
   {
   printf("alpha %f,delta %f, mag %d \n",stars[i].ra,stars[i].dec,stars[i].mag);
	}
  {
     short ref;
     long n;
     if (FSOpen("\psstgraphs.dat",0,&ref)) {
         printf("***ERREUR*** ouverture sstgraphs.dat\n");
         return;
     }
     n = sizeof(long);
     FSRead(ref, &n, (char*) &nposaxesimu);
   printf("nstarcatal = %d\n", nposaxesimu);
     n = nposaxesimu*sizeof(struct sstgraph);
     listord = (struct sstgraph*) malloc(n);
     FSRead(ref, &n, (char*) listord);
     FSClose(ref);
     printf("sstgraph.dat : read %d positions simulees\n", nposaxesimu);
  }
	printf("nombre d'etoiles selectionnees %d\n",listord[0].netdet);
	numaccro=0;flcatmin=1;
}
}

void exec_recons_sst()
{
int i,icour;
int cascour,casper,vudsper;
int valid;
float delbar,delebarre;
float deltan,deltintin;
float delcas,pocac;
sflxv=0.;
if(nfoundstars>0)
{
if(nfoundstars>MAXFOUNDSTARS) nfoundstars=MAXFOUNDSTARS;
laplushaute=0.;
for(i=0;i<nfoundstars;i++)
{
if(mfoundstars[i]>sflxv)
{sflxv=mfoundstars[i];
laplushaute=zfoundstars[i];
}
//printf("%d %f %f %f\n",i,tfoundstars[i],zfoundstars[i],mfoundstars[i]);
}

if(sflxv>CoupFlux)
{
maintenant=tfoundstars[0];
membar=laplushaute;
if(flagstar==0)
 {
 tdeb=maintenant;
 flagstar=1;
 }
 if(maintenant-tdeb>LowSeek)
 {
 vudsper=0;
 valid=0;
 delbar=0;
 delebarre=10000;
 deltan=0.;
 deltintin=1.;
for(i=0;i<nmemechant;i++)
{
if(sfluxvu[i]==0.) break;
separtps=maintenant-instant[i];
separenz=laplushaute-barette[i];
if(separtps>LowSeek && fabs(sfluxvu[i]-sflxv)/sflxv<.5)
 {
 if(separtps>tpsminfen && separtps<tpsmaxfen)
 {
 trace_point(fenetre_Find_Period);
 nptsfenper++;
 if(nptsfenper==QuorumRafrai && fenetre(fenetre_Find_Period))
 { selectgra(fenetre_Find_Period);efface(fenetre_Find_Period);nptsfenper=0;}
 }
 // recherche de periode
 if(separtps<HighSeek)
 {
 casper=(int)((Period-LowSeek)/StepCas);
 cascour=(int)((separtps-LowSeek)/StepCas);
 delcas=separtps-StepCas*(cascour+.5)-LowSeek;
 pocac=fabs(delcas/StepCas);
// printf("separtps %f delcas %f pocac %f\n",separtps,delcas,pocac);
 
 
 if(cascour>0 && cascour<ncasdiftps-1)
  {
  cases[cascour]+=25./separtps*(1.-pocac);
  if(delcas>0.){cases[cascour+1]+=25./separtps*pocac;}
  else {cases[cascour-1]+=25./separtps*pocac;}
  somcas+=25./separtps;
   if(fabs(cascour-casper)<2)
   {vudsper++;
   }
  if(cases[cascour]>maxcas)
   {
   maxcas=cases[cascour];
   poidsvois=cases[cascour-1]+cases[cascour]+cases[cascour+1];
   NewPeriod=LowSeek+
   (cascour+.5+(cases[cascour+1]-cases[cascour-1])/poidsvois)*StepCas;
    if(poidsvois>5*somcas/ncasdiftps){
    flagChPer=0;
    if(fabs(Period-NewPeriod)>1.)
    {nstastat=0;ndsper=0;tauxper=0.;
    { if(DonnePeriod() != -1 && fenetre(fenetre_Find_Period))
    {selectgra(fenetre_Find_Period);efface(fenetre_Find_Period);nptsfenper=0;}}}
    if(fabs(Period-NewPeriod)>.1) flagChPer=1;
    if(fabs(Period-NewPeriod)>.010)
    {raz_fen_phase();}
    Period=NewPeriod;
    }
//   printf("maxcas %f somcas/ncasdiftps %f Periode : %f\n",maxcas,
//   somcas/(float)(ncasdiftps),Period);
   }
  }
  }
  delbar=separenz;
  deltan=separtps-Period;
  if(DonnePeriod() != -1 
  && fabs(deltan-Period)<FourchTps && fabs(deltan-Period-2*deltintin)<FourchTps 
  && fabs(delebarre-2*delbar)<8 && valid==1)
  {
	tracephase();
//	membar=.3*barette[i]+.7*membar;
	nptsfenphs++;
	if(nptsfenphs==RafraiPhase)	{raz_fen_phase();}
//	printf("periode %8.3f phases tpsdiodes %8.1f sens %d\n",DonnePeriod(),DonnePhase(),DonneSens());

  }
  if(DonnePeriod() != -1 && fabs(deltan)<FourchTps && fabs(delbar)<8)
  {valid=1;delebarre=delbar;deltintin=deltan;}
  
 }
}
nstastat++;
if(vudsper!=0)ndsper++;
if(nstastat>StatValidePeriode)
{
	tauxper=(float)(ndsper)/(float)(nstastat);

//	printf("efficacite periode %f\n",tauxper);
}

if(nstastat==50)	{nstastat=0;ndsper=0;}
if(DonnePeriod()==-1){
raz_fen_phase();
if(fenetre(fenetre_Find_Period) && tpsminfen != 10.){
    				selectgra(fenetre_Find_Period);
    				graph->xmin=10.;
    				graph->xmax=40.;
    				graph->xpas=10.;
    				Liminffen=graph->xmin;
    				Limsupfen=graph->xmax;
				    tpsminfen=graph->xmin;
    				tpsmaxfen=graph->xmax;
    				passargu(Liminffen,Limsupfen,Papa);
    				efface(fenetre_Find_Period);
    				retrace(fenetre_Find_Period);
					nptsfenper=0;
}

} else
{
if(flagChPer==1){
		selectgra(fenetre_Find_Period);
    				graph->xmin=DonnePeriod()-1.;
    				graph->xmax=DonnePeriod()+1.;
    				graph->xpas=1.;
    				Liminffen=graph->xmin;
    				Limsupfen=graph->xmax;
				    tpsminfen=graph->xmin;
    				tpsmaxfen=graph->xmax;
    				passargu(Liminffen,Limsupfen,Papa);
    				efface(fenetre_Find_Period);
    				retrace(fenetre_Find_Period);
					nptsfenper=0;
}
}
if(fenetre(fenetre_Find_Period))
{selectgra(fenetre_Find_Period);ecritD(fenetre_Find_Period,FP_vperesti," %7.3f",DonnePeriod());}
if(Liminffen==LowSeek && DonnePeriod() != -1 && fenetre(fenetre_Find_Period)) 
{
    				selectgra(fenetre_Find_Period);
    				graph->xmin=DonnePeriod()-1.;
    				graph->xmax=DonnePeriod()+1.;
    				graph->xpas=1.;
    				Liminffen=graph->xmin;
    				Limsupfen=graph->xmax;
				    tpsminfen=graph->xmin;
    				tpsmaxfen=graph->xmax;
    				passargu(Liminffen,Limsupfen,Papa);
    				efface(fenetre_Find_Period);
    				retrace(fenetre_Find_Period);

}
}

icour=0;
pluvieux=maintenant;
for(i=nmemechant-1;i>0;i--)
 {
 icour++;
 if(sfluxvu[i-1]==0.) continue;
 if(instant[i]<pluvieux)pluvieux=instant[i];
 sfluxvu[i]=sfluxvu[i-1];
 instant[i]=instant[i-1];
 barette[i]=barette[i-1];
 phd0[i]=phd0[i-1];
  if(DonnePeriod()!=-1){if(maintenant-instant[i]>3.*DonnePeriod())
   {
    sfluxvu[i]=0.;
    instant[i]=0.;
    barette[i]=0.;
    phd0[i]=0.;
   }
   }
 }
sfluxvu[0]=sflxv;
instant[0]=maintenant;
barette[0]=laplushaute;
phd0[i]=0.;
{float peri=DonnePeriod();
if (peri>0.)phd0[i]=360./peri*(maintenant-peri*(int)(maintenant/peri));
}
// bloc recherche de l'axe ->identification des etoiles
{float difftps=0.;

if(DonnePeriod()>0. && (maintenant - pluvieux)>DonnePeriod())
{
nlistevue=0;
for(i=0;i<nmemechant-1;i++)
{
 if(sfluxvu[i]==0. |maintenant-instant[i]>Period) break;
difftps=(instant[i]-instant[i+1])*DPI/Period;
nlistevue++;
deltap[i]=difftps;
rapflx[i]=sfluxvu[i+1]/sfluxvu[i];
}
//printf("%d etoiles detectees par tour\n",nlistevue);
if(fenetre(fenetre_Find_Period))
{selectgra(fenetre_Find_Period);ecritD(fenetre_Find_Period,Find_net," %d",nlistevue);}
numaccro=1;
if(numaccro==0) chercheaccro(flcatmin);




}}
// fin du bloc recherche de l'axe ->identification des etoiles

}
}
}
void chercheaccro(int seuil)
{
/*
#define nstpturn 100
struct sstgraph {
  int netdet; // nombre d'etoiles detectees
  float raax; // coordonnees equatoriales de l'axe de rotation du telescope
  float decax;
  long numcat[nstpturn]; //numero de l'etoile dans le fichier de depart
  float phase[nstpturn]; //phase de l'etoile (angle/Nord - Est=+90ˇ - ds plan perp axe rot ballon)
  float cosinus[nstpturn]; // de l'angle entre la direction de l'etoile et l'axe de rotation du ballon
  float fluxcat[nstpturn]; //flux de l'etoile en unites arbitraires
  int	icl[nstpturn]; //ordre de l'etoile pour un classement par phases croissantes
};
*/
int i,j;int nstartour;float temps[nstpturn],fluxion[nstpturn];
nstartour=0;
for(i=0;i<nposaxesimu;i++)
	{
	for(j=0;j<listord[i].netdet-1;j++)
		{
		deltapcat[j]=listord[i].phase[j+1]-listord[i].phase[j];
		rapflxcat[j]=listord[i].fluxcat[j+1]/listord[i].fluxcat[j];
		}
	}


}
void trace_point(int fen)
{
if(fenetre(fen)) symbole(fen, separtps, sflxv, 2, rond, 0, bleu);
}
void tracephase(void)
{
float tempsplie,tempsmoy;float peri;float ecartps,ecarz;
peri=DonnePeriod();
if (peri>0.)
{
tempsplie=360./peri*(maintenant-peri*(int)(maintenant/peri));
if(!flagvuefiltree)
{
if(fenetre(fenetre_Phase_SST))
{symbole(fenetre_Phase_SST,tempsplie,separtps-2.*peri, 10, carre, 0, bleu);
symbole(fenetre_Phase_SST,tempsplie,.005*(separenz), 10, croix, 0, rouge);}
}

{
if(ncoups==0)memtemps=((int)(tempsplie)/pasenphase)*pasenphase;
 if(tempsplie-memtemps<pasenphase && tempsplie>=memtemps)
 {
  ncoups++;
  somtps+=tempsplie;
  somecartps+=separtps-2.*peri;
  somecarz+=(separenz)*.005;
 } else
 {
  if(ncoups>0)
  {
   tempsmoy=somtps/(float)(ncoups);
//   printf("temps,ecart : %f %f\n",tempsmoy,somecartps/ncoups);
   ecartps =somecartps/(float)(ncoups);
   ecarz = somecarz/(float)(ncoups);
   charge_calc_phase(tempsmoy,ecartps,ecarz);
if(fenetre(fenetre_Phase_SST) && flagvuefiltree)
{   symbole(fenetre_Phase_SST,tempsmoy,ecartps, 10, carre, 0, bleu);
   symbole(fenetre_Phase_SST,tempsmoy,ecarz, 10, croix, 0, rouge);
}
   ncoups=0;
   somtps=0;
   somecartps=0;
   somecarz=0;
   memtemps=((int)(tempsplie)/pasenphase)*pasenphase;
  {
  ncoups++;
  somtps+=tempsplie;
  somecartps+=separtps-2.*peri;
  somecarz+=(separenz)*.005;
 }
 }
 }

}
}}
float DonnePeriod(void)
{
if(tauxper>.4){
return Period;}
else return -1.;
}
int CpFl(void)
{
return (int)(CoupFlux);
}
void raz_fen_phase()
{
int i;
if (fenetre(fenetre_Phase_SST)){selectgra(fenetre_Phase_SST);
efface(fenetre_Phase_SST);}
memtemps=-1.;ncoups=0;nptsfenphs=0;phaseret=-9999.;Sens=-9999;
CasesVides=ncasphase;for(i=0;i<ncasphase;i++) statcalcphase[i]=0;

}
void charge_calc_phase(float tempsmoy,float ecartps,float ecarz)
{
int caserne;
caserne=(tempsmoy+pasenphase*.5)/pasenphase;
if(caserne==ncasphase)caserne=0;
if(caserne < ncasphase && caserne >= 0)
{
if(statcalcphase[caserne]<5 | (fabs(ecarmoytps[caserne]-ecartps)<.02 &&
fabs(ecarmoyz[caserne]-ecarz)<3.))
{ecarmoytps[caserne]*=statcalcphase[caserne];
ecarmoyz[caserne]*=statcalcphase[caserne];
ecarmoytps[caserne]+=ecartps;
ecarmoyz[caserne]+=ecarz;
statcalcphase[caserne]++;
ecarmoytps[caserne]/=statcalcphase[caserne];
ecarmoyz[caserne]/=statcalcphase[caserne];
CasesVides=calc_phase();
}}}
void tabulesincos(void)
{
int i;float pasenrad;
pasenrad=pasenphase*.01745329;
for(i=0;i<ncasphase;i++)
{
sini[i]=sin(pasenrad*(i));
cosini[i]=cos(pasenrad*(i));
}
}
int calc_phase(void)
{
int i,vides;float cocoz,sisiz,cocotps,sisitps;
float phasez,phasetps;float sensrot;
static float phasemem;
vides =0;
for(i=0;i<ncasphase;i++)
{
if(statcalcphase[i]==0) vides++;
}

cocoz=0;
sisiz=0;
cocotps=0;
sisitps=0;
for(i=0;i<ncasphase;i++)
{
if(statcalcphase[i]>0)
{
cocoz+=cosini[i]*ecarmoyz[i];
sisiz+=sini[i]*ecarmoyz[i];
cocotps+=cosini[i]*ecarmoytps[i];
sisitps+=sini[i]*ecarmoytps[i];
}
}
phasez=-atan2(cocoz,sisiz);
phasetps=-atan2(cocotps,sisitps);
sensrot=phasez-phasetps;
if(sensrot>3.14)sensrot-=6.28;
if(sensrot<-3.14)sensrot+=6.28;
if(sensrot>0.) Sens =1;
if(sensrot<=0.) Sens =-1;
//arbitraire, ajuste sur la simulation
Phase=phasez*57.29577;
phaseret=-9999.;
if(CasesVides<3 && fabs(Phase-phasemem)<10.){ phaseret=Phase;}
if(CasesVides<4) phasemem=Phase;
//printf("phases : numdiodes %f, tpsdiodes %f sens %d\n",phasez,phasetps,Sens);
return vides;}
float DonnePhase(void)
{
return phaseret;
}
int DonneSens(void)
{
if(CasesVides<3){return Sens;}else {return -9999;}
}