/*********************************************************************
 *                                                                   *
 * JOINSOURCES.C  -  Joins sources of same frequency                 *
 *                                                                   *
 * (c) 2004 by Bernd Vollmer, CDS, Observatoire de Strasbourg        *
 *                                                                   *
 *********************************************************************/
#include <stdio.h>
#include "nrutil.h" 
#include <math.h>
#include "code.h"

local int fcompare(bodyptr,bodyptr);

void joinsource(bodyptr p)
{
  int i,j, k, pl, ps, pb, strc;
  real sum, error, sum1, error1, factor;
  comppareptr q;
  bodyptr ppp[maxparents],sss[maxsons],bbb[maxbrothers];
  short pflag[maxparents],bflag[maxbrothers],sflag[maxsons];

  factor=1.2;                           /* fudge factor for flux error */
  for (i=0;i<Nbrothers(p);i++) bflag[i]=0;
  for (i=0;i<Nparents(p);i++) pflag[i]=0;
  for (i=0;i<Nsons(p);i++) sflag[i]=0;
  pl=ps=pb=0;

  if (Nparents(p) > 0) {                 /* search within Parents */
    for (i=0;i<Nparents(p);i++) ppp[i]=(bodyptr) Parents(p)[i];  /* remember parents */
    for (i=0;i<Nparents(p);i++) {
      if (Frequency(Parents(p)[i]) == Frequency(p)) {            /* at the same frequency */
	sum1=Flux(Parents(p)[i]);
	error1=Fluxerror(Parents(p)[i]);
	if (Flux(p)+factor*Fluxerror(p) > sum1-error1 && Flux(p)-factor*Fluxerror(p) < sum1+error1) {           /* if same flux source = brother */ 
	  pflag[i]=BROTHER;
	  pb++;
	}
	else if (sum1-error1 >= Flux(p)+factor*Fluxerror(p)) {  /* if greater flux source = parent */
	  pflag[i]=PARENT;
	  pl++;
	}
	else if (sum1+error1 <= Flux(p)-factor*Fluxerror(p)) {  /* if smaller flux skip source */
	  pflag[i]=0;
	  ps++;
	}
	else {
	  printf("problem in joinsources.c  1 \n");
	  getchar();
	}
      } else {
	pflag[i]=PARENT;    /* if not at the same frequency source = parent */
	pl++;
      }
    }
    if (Nbrothers(p)+pb > maxbrothers) {
      printf("too many brothers in joinsources.c \n");
      exit(0);
    }
    j=0;
    for (i=0;i<Nparents(p);i++) {      /* reconfigure Brothers */
      if (pflag[i] == BROTHER) {
	Brothers(p)[Nbrothers(p)]=Parents(p)[i];
	Nbrothers(p)++;
      } 
      else if (pflag[i] == PARENT) {   /* reconfigure Parents */
	Parents(p)[j]=(nodeptr) ppp[i];
	j++;
      }
    }
    Nparents(p)=j;
    if (pl != j) {
      printf("problem in joinsources.c 2 \n");
      getchar();
    }
  }
  
  pl=ps=pb=0;

  if (Nsons(p) > 0) {         /* search within sons */
    for (i=0;i<Nsons(p);i++) sss[i]=(bodyptr) Sons(p)[i];   /* remember sons */
    for (i=0;i<Nsons(p);i++) {
      if (Frequency(Sons(p)[i]) == Frequency(p)) {  /* at the same frequency */
	sum1=Flux(Sons(p)[i]);
	error1=Fluxerror(Sons(p)[i]);
	if (Flux(p)+factor*Fluxerror(p) > sum1-error1 && Flux(p)-factor*Fluxerror(p) < sum1+error1) {           /* if same flux source = brother */
	  sflag[i]=BROTHER;
	  pb++;
	}
	else if (sum1-error1 >= Flux(p)+factor*Fluxerror(p)) {   /* if greater flux skip source */
	  sflag[i]=0;
	  pl++;
	}
	else if (sum1+error1 <= Flux(p)-factor*Fluxerror(p)) {   /* if smaller flux source */
	  sflag[i]=SON;
	  ps++;
	}
	else {
	  printf("problem in joinsources.c 3 \n");
	  getchar();
	}
      } else {        /* if not at the same frequency source = son */
	sflag[i]=SON;
	ps++;
      }
    }
    if (Nbrothers(p)+pb > maxbrothers) {
      printf("too many brothers in joinsources.c \n");
      exit(0);
    }
    j=0;
    for (i=0;i<Nsons(p);i++) {      /* reconfigure Brothers */
      if (sflag[i] == BROTHER) {
	Brothers(p)[Nbrothers(p)]=Sons(p)[i];
	Nbrothers(p)++;
      } 
      else if (sflag[i] == SON) {   /* reconfigure Sons */
	Sons(p)[j]=(nodeptr) sss[i];
	j++;
      } 
    }
    if (ps != j) {
      printf("problem in joinsources.c 4 \n");
      getchar();
    }
    Nsons(p)=j;
  }
  

  pl=ps=pb=0;
	 
  if (Nbrothers(p) > 1) {        /* search within brothers */
    q=sorting;
    j=Nbrothers(p);
    for (i=0;i<Nbrothers(p);i++) {
      Compnum(q)=(bodyptr) Brothers(p)[i];
      Compfrequency(q)=Frequency(Brothers(p)[i]);
      Compi(q)=i;
      q++;
    }
    qsort(sorting,Nbrothers(p),sizeof(comppare),fcompare);   /* sort by frequency */
    q=sorting;
    for (i=0;i<Nbrothers(p);i++) {
      Brothers(p)[i] = (nodeptr) Compnum(q);
      q++;
    }
    for (i=0;i<Nbrothers(p);i++) {
      if (Frequency(Brothers(p)[i]) == Frequency(p)) {  /* at the same frequency */
	sum=sum1=Flux(Brothers(p)[i]);
	error=error1=factor*Fluxerror(Brothers(p)[i])*factor*Fluxerror(Brothers(p)[i]);
	k=0;
	for (j=i+1;j<Nbrothers(p);j++) {  /* add flux of all sources */
	  strc=strncmp(Name(Brothers(p)[i]),Name(Brothers(p)[j]),3);
	  if (strc == 0 && Frequency(Brothers(p)[j]) == Frequency(Brothers(p)[i])) {
	    sum+=Flux(Brothers(p)[j]);
	    error+=factor*Fluxerror(Brothers(p)[j])*factor*Fluxerror(Brothers(p)[j]);
	    k++;
	  }
	}
	error=sqrt(error);
	error1=sqrt(error1);
	if (Flux(p)+factor*Fluxerror(p) > sum1-error1 && Flux(p)-factor*Fluxerror(p) < sum1+error1) {     /* if same flux source = brother */
	  bflag[i]=BROTHER;
	} 
	if (Flux(p)+factor*Fluxerror(p) > sum-error && Flux(p)-factor*Fluxerror(p) < sum+error)  {        /* if sum flux = flux sources = brothers */
	  for (j=i+1;j<Nbrothers(p);j++) { 
	    strc=strncmp(Name(Brothers(p)[i]),Name(Brothers(p)[j]),3);
	    if (strc == 0 && Frequency(Brothers(p)[j]) == Frequency(Brothers(p)[i]))  {
	      bflag[j]=BROTHER;
	    }
	  }
	}
      } else bflag[i]=BROTHER;  /* not at the same frequency source = brother */
    }
  } else if (Nbrothers(p) == 1) {   /* case of single brother */
    sum1=Flux(Brothers(p)[0]);
    error1=factor*Fluxerror(Brothers(p)[0]);
    if (Frequency(Brothers(p)[0]) == Frequency(p) && Flux(p)+factor*Fluxerror(p) > sum1-error1 && Flux(p)-factor*Fluxerror(p) < sum1+error1) {
      bflag[0]=BROTHER;
    }
    if (Frequency(Brothers(p)[0]) != Frequency(p)) bflag[0]=BROTHER;
  }
  j=0;
  for (i=0;i<Nbrothers(p);i++) bbb[i]=(bodyptr) Brothers(p)[i];

  for (i=0;i<Nbrothers(p);i++) {    /* reconfigure Brothers */
    if (bflag[i] == BROTHER) {
      Brothers(p)[j]=(nodeptr) bbb[i];
      j++;
    }
  }
  
  if (j > 0) njoin++;
  Nbrothers(p)=j;
}

/*
 *  quicksort subroutine
 */ 
	
	 
local int fcompare(bodyptr a,bodyptr b)
{
  if ((Compfrequency(a)-Compfrequency(b)) > 0.0) return(1);
  else if ((Compfrequency(a)-Compfrequency(b)) < 0.0) return(-1);
  else return(0);
}