/* THIS PROGRAM IDENTIFIES THE MEMBERS OF A FAMILY WHICH ARE ASSOCIATED TO
   A BODY *idesig_to_start* AT A LEVEL *limit* (TO BE SET BELOW) 
   WORKS ON ANALYTICAL ELEMENTS OF MILANI AND KNEZEVIC, HCM OF ZAPPALA ET AL. */

#include<stdio.h>
#include<math.h>
#include<stdlib.h>

#include "nrutil.h"
#include "nrutil.c"

#define PI2 6.28318530717959 
#define AU 1.49597870e11 // IN METERS

#define GAUSS 0.01720209895
#define GM (GAUSS*GAUSS*365.25*365.25) //UNITS: JULIAN YEARS, AU, SUN MASS

#define FACTOR (AU/365.25/86400.) // CONVERSION FACTOR FROM AU/YR TO M/S

#define MAXN_IN_LIMIT 500 // MAXIMUM NUMBER OF BODIES ALLOWED WITHIN A *LIMIT* ABOUT ANY BODY 
                           // IF EXCEEDED, PROGRAM GIVES A RUN-TIME ERROR MESSAGE AND EXISTS  

static double sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)

void fgetfield(FILE *fp,char s[],int length);

int main(void){

  char field[100];
  int RFL,QCM,QCO,imy;
  long int i,j,k,l,inew,ide,ntpnumb,ntpmult,lastnum,ntp;
  long int *idesig,*nassoc,**assoc,nfamily,*family;
  long int idesig_to_start, k_to_start;
  double *sema,*ecc,*sinc,n,g,s,*mag,lyap;
  double limit,sqrfac,metrics,maglimit;
  double elow,ehigh;
  FILE *fp;

  //********************************************************************************

  scanf("%ld %lg",&idesig_to_start,&limit);

  //  printf("Enter id of the central body (e.g. 490 is Veritas) : ");
  // scanf("%ld",&idesig_to_start);
  //  idesig_to_start = 158; // ENTER HERE THE DESIGNATION NUMBER OF THE BODY YOU WISH 
                         // TO ASSOCIATE THE FAMILY WITH 
  //  printf("Enter cut-off : ");
  // scanf("%lg",&limit);
  //  printf("Enter limiting magnitude of bodies to be considered: ");
  //  scanf("%lg",&maglimit);
  maglimit = 100;
  // limit = 100; // SEPARATION LIMIT, IN M/S, MUST BE LARGER THAN QRL (ZAPPALA ET AL.)
  
  limit = limit*limit;

  //********************************************************************************

  sqrfac = SQR(FACTOR); // TO AVOID SQUARE ROOT IN METRICS, WORK WITH METRICS^2  

  ntpnumb = 302212;
  lastnum = 312935;
  ntpmult = 0;

  ntp = ntpnumb + ntpmult;
  
  sema = dvector(1,ntp); ecc = dvector(1,ntp); sinc = dvector(1,ntp);
  mag =  dvector(1,ntp);
  idesig = lvector(1,ntp);

  nassoc = lvector(1,ntp);
  assoc = lmatrix(1,ntp,1,MAXN_IN_LIMIT);
  family = lvector(1,ntp);

  // READ THE PROPER ELEMENTS OF NUMBERED ASTEROIDS 
  fp = fopen("numb.syn","r");

  while(fgetc(fp) != '\n'); 
  while(fgetc(fp) != '\n');

  for(i=1;i<=ntpnumb;i++){
      
    fgetfield(fp,field,10);
    idesig[i] = atol(field);

    fscanf(fp,"%lg",&mag[i]);    
    fscanf(fp,"%lg %lg %lg %lg %lg %lg",
       &sema[i],&ecc[i],&sinc[i],&n,&g,&s);
    fscanf(fp,"%lg %d",&lyap,&imy);

    //    printf("%ld %lg %lg %lg\n",idesig[i],sema[i],ecc[i],sinc[i]);
  }
  fclose(fp);


  // WHAT IS THE BODY TO START WITH ? 
  i = 1;
  while(idesig[i++] != idesig_to_start && i<=ntp);

  if(i>ntp){ 
    fprintf(stderr,"The central body, %ld, not found ...", idesig_to_start);
    fprintf(stderr,"Exiting now ...");
    return 0;
  }

  k_to_start = i-1;

  // NASSOC[I] IS # OF BODIES WITHIN *LIMIT* DISTANCE TO BODY *I*, INITIALLY SET TO ZERO
  for(i=1;i<=ntp;i++)
    nassoc[i] = 0; 
  
  // SELECTION OF BODIES LINKED WITH *idesig_to_start*
  k = k_to_start;         
  nfamily = 1;
  family[nfamily] = k;       // add *idesig_to_start* as 1st body to the family list 

  i = 1;  /* run with *i* over the current list of family members, 
       attaching additional branches */
  
  while(i <= nfamily){  // stop if *i* is larger than the current family list: no more branches
   
    k = family[i];   /* run over the bodies associated to *i* body */
          
    for(j=1;j<=ntp;j++)
      
      if(j != k){
  
       // D1 METRICS SQUARED (!) OF ZAPPALA ET AL. 1994
       metrics = 
         1.25*SQR( 2.00*(sema[j] - sema[k])/(sema[k] + sema[j]) ) + 
         2.00*SQR(  ecc[j]  - ecc[k]           ) +
         2.00*SQR(  sinc[j] - sinc[k]          );    
  
       metrics = 2. * metrics*GM/(sema[j] + sema[k]) * sqrfac;
  
       // IF J-TH BODY IS CLOSER THAN *LIMIT*, ASSOCIATE IT TO I-TH BODY

       // diagonal cuts for Flora
         elow = 0.786720*sinc[j] + 0.0438681;
         ehigh =  1.17162*sinc[j] + 0.0658839; 

         if(metrics <= limit && mag[j] <= maglimit){
       // if(metrics <= limit && mag[j] <= maglimit && ecc[j]>elow && ecc[j]<ehigh){ //diagonal cut for Flora
       //if(metrics <= limit && mag[j] <= maglimit && sinc[j]<0.031759651){ // cut for Massalia
       //if(metrics <= limit && mag[j] <= maglimit && sinc[j]>0.031759651){ // cut for Nysa-Polana
         nassoc[k] ++;
    
           if(nassoc[k] > MAXN_IN_LIMIT){
           fprintf(stderr,"Number of associated bodies to %ld exceeds MAXN_IN_LIMIT ...\n",k);
           fprintf(stderr,"Increase MAXN_IN_LIMIT, Exiting ...");
             return 0;
           }
    
           assoc[k][nassoc[k]] = j;
           } 

      }      
        
    for(j=1;j<=nassoc[k];j++){   
      
      inew = 1;      /* see if bodies associated to *k* are already in the family list */
      for(l=1;l<=nfamily;l++)
      if(assoc[k][j] == family[l] )
         inew = 0;

      if(inew == 1){    /* if this is a new body, add it to the family list */
       nfamily++;
       fprintf(stderr,"%ld\n",nfamily);
       family[nfamily] = assoc[k][j];
      }
  
    }

    i++;
  }

  // OUTPUT THE FAMILY LIST 
  fp = fopen("family.list","w");
  for(i=1;i<=nfamily;i++){

    ide = family[i];
    //    fprintf(fp,"%6ld %6ld %8lg %7lg %7lg %6lg %s\n",
    //       i,idesig[ide],sema[ide],ecc[ide],sinc[ide],mag[ide],sdesig[ide]); 
    fprintf(fp,"%6ld %6ld %8lg %7lg %7lg %6lg\n",
            i,idesig[ide],sema[ide],ecc[ide],sinc[ide],mag[ide]); 
  }
  fclose(fp);

  return 1;
}

void fgetfield(FILE *fp,char s[],int length){
  int i;
  
  for(i=0;i<length;i++)    
    s[i]=fgetc(fp);
  s[i]='\0';
}