/*****************************************************************
 *                                                               *
 *  SPECFIND main routine                                        *
 *                                                               *
 *  (c) 2004 by Bernd Vollmer, CDS, Observatoire de Strasbourg   *
 *                                                               *
 *****************************************************************/
#define global
#include <stdio.h>
#include "nrutil.h"
#include "code.h"



long zeit11,zeit22; 
real em1,em2;
int totalnum1, totalnum2, totalnum3, totalnum4, totalnum5, nnbody, jj, totalspread;

void main()
{
	void odeint(float[],int,float,float,float,float,float,int *,int *,void (*derivs)(float,float[],float[]),
		void(*bsstep)(float[],float[],int,float *,float,float,float[],float *,float *,
		void(*)(float,float[],float[])));
	FILE *datei_ptr;
	int iii,jjj;
	char Datei[13];
	float *ystart;
	float h1,ttiny=1.e-3;
	float x1=0.0,x2=8.0e8,hmin=1.0e-5,epso=0.1;
	
	int counter,i, j, newdim, totalnum, filenum, nspectra1;
	real distance,zwischen,bilanz,dist,dist0;
	real am1,am2,am3,em3,hhelp1,hhelp2,hhelp3,hhelp4,hhelp5,hhelp6;
	long zeit1,zeit2,zeit3;
	bool ask0;
	vector rrr,sss,ttt,xa1,xa2,xa3,va1,va2,va3,xe1,xe2,xe3,ve1,ve2,ve3;
	bodyptr p,pq;
	catalogueptr rr;
	char cc[10],dd[19],ee[10],aaa[1000],bbb[1000],ccc[1000];

	options="bh86";                     /* Barnes&Hut criterion for tree */
	CLRV(offset);	                    /* tree offset */
	theta=0.0;

	maxbrothers=50;                     /* maximum number of brothers */
	maxparents=50;                      /* maximum number of parents  */
	hhelp1=MAX(maxbrothers,maxparents);
	maxsons=50;                         /* maximum number of brothers */
	hhelp1=MAX(hhelp1,maxsons);          
      	
      
	time(&zeit11);

	maketr=0;                     /* for Brothers in maketree */
	nspectra1=0;                  /* number of spectra found */
	rsize=1.e3;                  /* initial root size for tree */
	nbody=0;                     /* total number of sources */
	nnbody=0;

	cats=(catalogueptr) calloc(1000,sizeof(catalog));   /* maximum number of catalogues */

	/* begin reading of catalogue file */
	  
	totalnum=totalnum1=totalnum2=totalnum3=totalnum4=totalnum5=0;
	datei_ptr=fopen("catalogues.dat","r");
	fscanf(datei_ptr,"%d \n",&iii);	
	printf("number of catalogues: %i \n",iii); 
	ncatalogue=iii;
	maxstr=0.;
	for (rr=cats;rr<cats+ncatalogue;rr++) {
	  fscanf(datei_ptr,"%s %s %f %f %f %f %s \n",&cc,&ee,&am1,&am2,&am3,&em1,&dd);
	  sprintf(Cname(rr),cc);
	  sprintf(Cbib(rr),dd);
	  sprintf(Cacro(rr),ee);
	  Cfrequency(rr)=am1;
	  Cresolution(rr)=am2;
	  Cmaxstructure(rr)=am3;
	  Csnumber(rr)=em1;
	  totalnum+=em1;
	  maxstr=FMAX(maxstr,am3);
	}	
	fclose(datei_ptr);

	/* end reading of catalogue file */

	maxstr=maxstr*60.;
	filenum=totalnum/3.e5+0.9999;  /* calculate number of sphere subdivisions */
	printf("number of output files %i \n",filenum);

	dist0=1.296e6;
	for (j=1;j<=filenum+2;j++) {    /* loop over subdivisions */
	  
	  if (j <= filenum) pole=0;
	  else if (j == filenum+1) pole=1;
	  else if (j == filenum+2) pole=2;
	  
	  /*if (j == 1) pole=1;
	    else if (j == 2) pole=2;
	    else getchar();*/
	  
	  if (pole > 0) printf("\n pole search \n");
	  
	  nbody=0;
	  
	  bodytab=(bodyptr) calloc(600000,sizeof(body));  /* maximum number of sources */
	  sorting=(comppareptr) calloc(1000,sizeof(comppare));   /* maximum number of points in spectrum */
	  
	  am1=em1=(j-1)*dist0/filenum;
	  am2=em2=j*dist0/filenum;
	  if (am2 > dist0) am2=em2=dist0;
	  minborder=am1-3.*maxstr;
	  maxborder=am2+3.*maxstr;
	  if (pole > 0) {
	    am1=em1=minborder=0.0;
	    am2=em2=maxborder=dist0;
	  }
	  printf("\n \n");
	  printf("RA from %5.3e to %5.3e arsec \n",minborder,maxborder);
	  printf("\n");
	  
	  /* begin reading of catalogues */
	  
	  for (rr=cats;rr<cats+ncatalogue;rr++) {
	    sprintf(dateiname,"%sout.dat",Cname(rr));
	    printf("reading catalogue %s \n",Cname(rr));
	    resol=Cresolution(rr);
	    config();  
	    Csnumber(rr)=zaehler;
	  }
	  
	  /* end reading of catalogues */ 

	  pnumber=bnumber=snumber=0;
	  
	  printf("begin maketree with %i particles \n",nbody);
	  excor=cmore=0;
	  maketree(bodytab,nbody);       /* built up of the tree */
	  
	  printf("maketree done...\n");
	  
	  for (p=bodytab;p<bodytab+nbody;p++) hackgrav(p, Mass(p) > 0.0);    /* proximity search */	  

	  /* begin check for RA = 00 00 00 */

	  doedge=FALSE;
	  if (j == 1 || j == filenum) {
	    doedge=TRUE;
	    dist0=1.296e6;
	    for (p=bodytab;p<bodytab+nbody;p++) {
	      if (Pos(p)[0] > dist0-3.*maxstr) {
		Mass(p)=2.;
		Pos(p)[0]-=dist0;
	      } 
	      else if (Pos(p)[0] < 3.*maxstr) {
		Mass(p)=3.;
		Pos(p)[0]+=dist0;
	      }
	      else {
		Mass(p)=0.0;
	      }
	    }
	    maketree(bodytab,nbody);       /* built up of the tree */
	    printf("maketree done...\n");
	    for (p=bodytab;p<bodytab+nbody;p++) {
	      if (Mass(p) == 2.) Pos(p)[0]+=dist0;
	      if (Mass(p) == 3.) Pos(p)[0]-=dist0;
	      if (Mass(p) > 0.0) hackgrav(p, Mass(p) > 0.0);    /* proximity search */
	      if (Mass(p) == 2.) Pos(p)[0]-=dist0;
	      if (Mass(p) == 3.) Pos(p)[0]+=dist0;
	    }
	    for (p=bodytab;p<bodytab+nbody;p++) {
	      if (Mass(p) == 2.) Pos(p)[0]+=dist0;
	      if (Mass(p) == 3.) Pos(p)[0]-=dist0;
	      Mass(p)=1.0;
	    }
	  }  
	  
	  /* end check for RA = 00 00 00 */

	  time(&zeit22);
	  printf("time used for search: %ld sec \n",zeit22-zeit11);
	  
	  printf(" %i parents  %i brothers  %i sons found\n",pnumber,bnumber,snumber);
	  printf(" %i multiple coordinates found \n",excor);


	  /* begin of source joining */
	  
	  time(&zeit11);
	  njoin=0;
	  for (p=bodytab;p<bodytab+nbody;p++) {
	    joinsource(p);
	  }
	  time(&zeit22);
	  printf(" %i brother sources joined \n",njoin);
	  printf("time used for joining: %ld sec \n",zeit22-zeit11);
	  
	  /* end of source joining */

	  /* begin of source crosscheck */
	  
	  time(&zeit11);
	  ncross=nresolved=0;
	  for (p=bodytab;p<bodytab+nbody;p++) {
	    crosssource(p);
	  }
	  time(&zeit22);
	  printf("%i sources dropped when cross checked \n",ncross);
	  printf("time used for crossing: %ld sec \n",zeit22-zeit11);
	  printf("%i resolved sources \n",nresolved);
	  
	  /* end of source crosscheck */

	  /* begin of source compression */
	  
	  time(&zeit11);
	  ncompress=ndouble=0;
	  for (p=bodytab;p<bodytab+nbody;p++) {
	    compresssource(p);
	  }
	  time(&zeit22);
	  printf("%i sources dropped when compressed \n",ncompress);
	  printf("%i double sources found \n",ndouble);
	  printf("time used for compressing: %ld sec \n",zeit22-zeit11);
	  
	  /* end of source compression */

	  /* Change flux error */
	  for (p=bodytab;p<bodytab+nbody;p++) if (Fluxerror(p) < 0.2*Flux(p)) Fluxerror(p)=0.2*Flux(p);
	  
	  /* begin of spectrum extraction */
	  
	  time(&zeit11);
	  nspectra=0;
	  for (p=bodytab;p<bodytab+nbody;p++) {
	    specrun=0;
	    spectra(p);
	  }
	  time(&zeit22);
	  printf("time used for spectra: %ld sec \n",zeit22-zeit11);
	  printf("number of spectra found: %i \n",nspectra);
	  nspectra1+=nspectra;

	  /* end of spectrum extraction */

	  /* begin of spectra check */ 

	  speccheck();
	  printf("check of spectra done \n ");

	  /* end of spectra check */


	  /* begin of spectra spreading */

	  spreadspec();
	  printf("number of spread spectra: %i \n",totalspread);

	  /* end of spectra spreading */	   

	  /* begin of output */
	  
	  jj=j;
	  output();

	  /* end of output */

	  free(bodytab);
	  free(sorting);
	  
	}

	free(cats);

	printf("\n \n");
	printf("total number of sources: %i \n",totalnum);
	printf("processed sources: %i \n",nnbody);
	printf("number of output sources: %i \n",totalnum1);
	printf("total number of spectra found: %i \n",nspectra1);
	printf("final number of spectra found: %i \n",totalnum2);
	printf("total number of associations: %i \n",totalnum3);
	printf("total number of independent spectra: %i \n",totalnum4);
	printf("total number of independent flat spectra: %i \n",totalnum5);
	printf("\n \n");

}