/*********************************************************************
 *                                                                   *
 * CROSSSOURCES.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"



void crosssource(bodyptr p)
{
  int i,j, k, pl, ps, pb, strc, str1, str2;
  bodyptr ppp[maxparents],sss[maxsons],bbb[maxbrothers];
  short pflag[maxparents],bflag[maxbrothers],sflag[maxsons];
  real sum, error, factor=0.8;

  str1=0;
  if (Nbrothers(p) > 0) {                        /* checks brothers */
    for (i=0;i<Nbrothers(p);i++) bbb[i]=(bodyptr) Brothers(p)[i];   /* remember brothers */
    for (i=0;i<Nbrothers(p);i++) {
      strc=strcmp(Name(p),Name(Brothers(p)[i])); /* does name already exist? */
      if (strc == 0) {                           /* if yes skip */ 
	bflag[i]=0;
	str1++;
	ncross++;
      } else bflag[i]=BROTHER;                  /* if no keep as brother */
    }
    if (str1 > 0) {
      bflag[i]=0;
    }
    j=0;
    for (i=0;i<Nbrothers(p);i++) {              /* reconfigure Brothers */
      if (bflag[i] == BROTHER) {
	Brothers(p)[j]=(nodeptr) bbb[i];
	j++;
      }
    }
    Nbrothers(p)=j;
  }
  
  if (Nparents(p) > 0) {                          /* check 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(p) == Frequency(Parents(p)[i])) {  /* at the same frequency */
	if (Nsons(Parents(p)[i]) == 0) pflag[i]=0;
	else if (Nsons(Parents(p)[i]) > 0) {      /* source = son of parent? */
	  str1=0;
	  for (j=0;j<Nsons(Parents(p)[i]);j++) {
	    strc=strcmp(Name(p),Name(Sons(Parents(p)[i])[j]));
	    if (strc == 0) str1+=1;
	  }
	  if (str1 >= 1) {
	    pflag[i]=PARENT;
	  } else if (str1 > 1) {
	    printf("problem with parents in crosssources.c \n");
	    getchar();
	  }	
	  else {
	    pflag[i]=0;
	    ncross++;
	  }
	}
      } else pflag[i]=PARENT;  /* not at the same frequency source = parent */
    }
    j=0;
    for (i=0;i<Nparents(p);i++) {                  /* reconfigure Parents */
      if (pflag[i] == PARENT) {
	Parents(p)[j]=(nodeptr) ppp[i];
	j++;
      }
    }
    Nparents(p)=j;
  }

  
  if (Nsons(p) > 0) {                              /* check 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(p) == Frequency(Sons(p)[i])) { /* at the same frequency */
	if (Nparents(Sons(p)[i]) == 0) sflag[i]=0;
	else if (Nparents(Sons(p)[i]) > 0) {
	  str1=0;
	  for (j=0;j<Nparents(Sons(p)[i]);j++) {  /* source = parent of son? */
	    strc=strcmp(Name(p),Name(Parents(Sons(p)[i])[j]));
	    if (strc == 0) str1+=1;
	  }
	  if (str1 >= 1) {
	    sflag[i]=SON;
	  } else {
	    sflag[i]=0;
	    ncross++;
	  }
	}
      } else sflag[i]=SON;  /* not at the same frequency source = son */
    }
    j=0;
    for (i=0;i<Nsons(p);i++) {                     /* reconfigure Sons */
      if (sflag[i] == SON) {
	Sons(p)[j]=(nodeptr) sss[i];
	j++;
      }
    }
    Nsons(p)=j;
  }


  if (Nbrothers(p) > 0) {                               /* check brothers */
    for (i=0;i<Nbrothers(p);i++) bbb[i]=(bodyptr) Brothers(p)[i];  /* remember brothers */
    for (i=0;i<Nbrothers(p);i++) {
      if (Frequency(p) == Frequency(Brothers(p)[i])) {  /* at the same frequency */
	if (Nbrothers(Brothers(p)[i]) == 0) bflag[i]=0;
	else if (Nbrothers(Brothers(p)[i]) > 0) {
	  str1=0;
	  for (j=0;j<Nbrothers(Brothers(p)[i]);j++) {   /* source = brother of brother? */
	    strc=strcmp(Name(p),Name(Brothers(Brothers(p)[i])[j]));
	    if (strc == 0) str1+=1;
	  }
	  if (str1 != 0) {
	    bflag[i]=BROTHER;
	  } 
	  else {
	    bflag[i]=0;
	    ncross++;
	  }
	}
      } else bflag[i]=BROTHER; /* not at the same frequency source = brother */
    }
    j=0;
    for (i=0;i<Nbrothers(p);i++) {                   /* reconfigure Brothers */
      if (bflag[i] == BROTHER) {
	Brothers(p)[j]=(nodeptr) bbb[i];
	j++;
      }
    }
    Nbrothers(p)=j;
  }
  
  /* begin check if source is resolved */

  if (Nsons(p) > 1) {
    for (i=0;i<Nsons(p);i++) if (Frequency(p) == Frequency(Sons(p)[i])) sflag[i]=0;
    else sflag[i]=-1;
    for (i=0;i<Nsons(p);i++) {
      if (sflag[i] == 0) {
	sflag[i]=1;
	pl=1;
	sum=Flux(Sons(p)[i]);
	error=Fluxerror(Sons(p)[i])*Fluxerror(Sons(p)[i]);
	for (j=i+1;j<Nsons(p);j++) {
	  if (Frequency(p) == Frequency(Sons(p)[j])) {
	    strc=strncmp(Name(Sons(p)[i]),Name(Sons(p)[j]),3);
	    if (strc == 0) {
	      sflag[j]=1;
	      sum+=Flux(Sons(p)[j]);
	      error+=Fluxerror(Sons(p)[j])*Fluxerror(Sons(p)[j]);
	      pl++;
	    }
	  }
	}
	error=sqrt(error);
	if (sum+error >= factor*Flux(p)-Fluxerror(p) && sum-error <= Flux(p)+Fluxerror(p)&& pl > 1) {
	  Wflag(p)=2;
	  nresolved++;
	}
      }
    }
  }
  
  /* end check if source is resolved */

}