/*****************************************************************
 *                                                               *
 * SPREADSPEC: Spreads spectra.                                  *
 *                                                               *
 * (c) 2004 by Bernd Vollmer, CDS, Observatoire de Strasbourg    *
 *                                                               *
 *****************************************************************/
#include "code.h"

local int diist;

void twofreq(bodyptr p) 
{
  int jj, i, j, kk, ninit, flag[50], ii;
  real freqhelp[51];
  bodyptr init[51];

  for (i=0;i<Nbrothers(p);i++) init[i]=(bodyptr) Brothers(p)[i];
  jj=1;
  freqhelp[0]=Frequency(p);
  for (i=0;i<Nbrothers(p);i++) {
    j=0;
    for (kk=0;kk<jj;kk++) if (Frequency(Brothers(p)[i]) == freqhelp[kk]) j=1;
    if (j == 0) {
      freqhelp[jj]=Frequency(Brothers(p)[i]);
      jj++;
    }
  }
  ninit=jj;
  
  if (ninit <= 2) {
    for (i=0;i<Nbrothers(p);i++) if (Frequency(p) == Frequency(Brothers(p)[i])) flag[i]=1;
    else flag[i]=0;
    j=0;
    for (i=0;i<Nbrothers(p);i++) {
      if (flag[i] == 1) {
	Brothers(p)[j]=(nodeptr) init[i];
	j++;
      }
    }
    Nbrothers(p)=j;
    Specindx(p)=0.0;
  } 
}


void distance(bodyptr pr, bodyptr pp)
{
  int i;
  real helpp, helpq, helph;
  bodyptr pz;

  diist=0;
  if (Majaxis(pr) == 0) helpp=Resolution(pr);
  else helpp=Majaxis(pr);
  if (Majaxis(pp) == 0) helpq=Resolution(pp);
  else helpq=Majaxis(pp);
  helpp+=helpq;
  helpp/=2.;
  DISTV1(helph,Pos(pr),Pos(pp));
  if (helph > helpp) diist=1;
  for (i=0;i< Nbrothers(pp);i++) {
    pz=(bodyptr) Brothers(pp)[i];
    if (Majaxis(pr) == 0) helpp=Resolution(pr);
    else helpp=Majaxis(pr);
    if (Majaxis(pz) == 0) helpq=Resolution(pz);
    else helpq=Majaxis(pz);
    helpp+=helpq;
    helpp/=2.;
    DISTV1(helph,Pos(pr),Pos(pz));
    if (helph > helpp) diist=1;
  }
}

extern int totalspread;

void spreadspec(void)
{	
  bodyptr p;
  int i, ij, ii, j, k, jj, kk, ninit, nfinal[51], nn, nn1, nn2, nn3, strc, strc1, qwe, flag[50], nnfinal, doing, pll, iz, npoints, in, zaq;
  bodyptr init[51], final[50], pq, pqq, pp, pr, pz;
  real freqhelp[51], freqhelp1[51], numin1, numin2, numax1, numax2, helpp, estflux[51], help, help1, dx, helph, helpq, resold;
  char aaa[12];

  for (p=bodytab;p<bodytab+nbody;p++) twofreq(p); /* security: remove source with only two different frequencies */
  
  for (p=bodytab;p<bodytab+nbody;p++) { /* find number of different frequencies for source p */
    for (i=0;i<Nbrothers(p);i++) init[i]=(bodyptr) Brothers(p)[i];
    jj=1;
    freqhelp[0]=Frequency(p);
    for (i=0;i<Nbrothers(p);i++) {
      j=0;
      for (kk=0;kk<jj;kk++) if (Frequency(Brothers(p)[i]) == freqhelp[kk]) j=1;
      if (j == 0) {
	freqhelp[jj]=Frequency(Brothers(p)[i]);
	jj++;
      }
    }
    ninit=jj;
    nn1=jj;
    nn2=-1;
    help=1.0e21;
    nn3=-1;
    pz=p;

    for (ij=0;ij<Nbrothers(p);ij++) { /* find maximum number of sources with different frequencies of all brothers */
      pq=(bodyptr) Brothers(p)[ij];
      j=0;
      for (i=0;i<Nbrothers(pq);i++) {
	strc=strcmp(Name(p),Name(Brothers(pq)[i]));
	if (strc == 0) j=1;
      }
      if (j == 1) {
	jj=1;
	freqhelp1[0]=Frequency(pq);
	for (i=0;i<Nbrothers(pq);i++) {
	  j=0;
	  for (kk=0;kk<jj;kk++) if (Frequency(Brothers(pq)[i]) == freqhelp1[kk]) j=1;
	  if (j == 0) {
	    freqhelp1[jj]=Frequency(Brothers(pq)[i]);
	    jj++;
	  }
	}    
	nfinal[ij]=jj;
	if (jj > nn1) {
	  nn1=jj;
	  nn2=ij;
	  pz=(bodyptr) pq;
	}
	helpp=Absz(pq)+Specindx(pq)*log10(Frequency(pq));
	estflux[ij]=fabs(Flux(pq)-pow(10.,helpp));
	if (estflux[ij] < help) {
	  help=estflux[ij];
	  nn3=ij;
	  pz=(bodyptr) pq;
	}
      }
    }
    
    if (nn2 != nn3) nn2=nn3;
    
    helpp=Absz(pz)+Specindx(pz)*log10(Frequency(p));
    help=fabs(Flux(p)-pow(10.,helpp));
    if (Specindx(p) == 0.0 && nn2 >= 0) { /* if a brother has a larger number of different frequencies then replace */
      pll=0;
      pp=pq;
      if (Majaxis(p) == 0) helpp=Resolution(p);
      else helpp=Majaxis(p);
      if (Majaxis(pp) == 0) helpq=Resolution(pp);
      else helpq=Majaxis(pp);
      helpp+=helpq;
      helpp/=2.;
      DISTV1(helph,Pos(p),Pos(pp));
      if (helph > helpp) pll=1;
      for (i=0;i< Nbrothers(pq);i++) {
	pp=(bodyptr) Brothers(pq)[i];
	if (Majaxis(p) == 0) helpp=Resolution(p);
	else helpp=Majaxis(p);
	if (Majaxis(pp) == 0) helpq=Resolution(pp);
	else helpq=Majaxis(pp);
	helpp+=helpq;
	helpp/=2.;
	DISTV1(helph,Pos(p),Pos(pp));
	if (helph > helpp) pll=1;
      }
      if (pll == 0) { 
	totalspread++;
	j=k=0;
	for (k=0;k<Nbrothers(p);k++) if (Frequency(Brothers(p)[k]) == Frequency(p)) j++;
	pq=pqq=(bodyptr) Brothers(p)[nn2];
	jj=0;
	for (i=0;i<Nbrothers(pq);i++) {
	  strc=strcmp(Name(p),Name(Brothers(pq)[i]));
	  if (strc == 0) jj=1;
	}
	if (jj == 1) {
	  if (Frequency(pq) != Frequency(p)) {
	    DISTV1(helph,Pos(p),Pos(pq));
	    if (helph <= (Resolution(p)+Resolution(pq))/2.) {
	      Brothers(p)[j]=(nodeptr) pq;
	      j++;
	      k++;
	    }
	  }
	  for (ii=0;ii<Nbrothers(pq);ii++) {
	    if (Frequency(p) != Frequency(Brothers(pq)[ii])) {
	      DISTV1(helph,Pos(p),Pos(Brothers(pq)[ii]));
	      if (helph <= (Resolution(p)+Resolution(Brothers(pq)[ii])/2.)) {
		Brothers(p)[j]=Brothers(pq)[ii];
		j++;
		k++;
	      }
	    }
	  }
	  qwe=k;
	  Nbrothers(p)=j;
	  Specindx(p)=Specindx(pq);
	  Absz(p)=Absz(pq);
	  Resolved(p)=Resolved(pq);
	}	
      }
    }  
  }
  
  for (p=bodytab;p<bodytab+nbody;p++) {  /* check for brothers of brothers */
    if (Nbrothers(p) > 0) {
      for (ii=0;ii<Nbrothers(p);ii++) flag[ii]=0;
      for (ii=0;ii<Nbrothers(p);ii++) {
	pq=(bodyptr) Brothers(p)[ii];
	if (Nbrothers(pq) > 0) {
	  if (fabs(Specindx(p)-Specindx(pq)) < 0.3) flag[ii]=1;
	}
      }
      for (ii=0;ii<Nbrothers(p);ii++) if (flag[ii] == 1) {
	pq=(bodyptr) Brothers(p)[ii];
	for (k=0;k<Nbrothers(pq);k++) {
	  nn=0;
	  strc=strcmp(Name(p),Name(Brothers(pq)[k]));
	  if (strc == 0) nn=1; 
	  for (ij=0;ij<Nbrothers(p);ij++) {
	    strc=strcmp(Name(Brothers(p)[ij]),Name(Brothers(pq)[k]));
	    if (strc == 0) nn=1;
	  }
	  if (nn == 0) {
	    pll=0;
	    pp=(bodyptr) Brothers(pq)[k];
	    if (Majaxis(p) == 0) helpp=Resolution(p);
	    else helpp=Majaxis(p);
	    if (Majaxis(pp) == 0) helpq=Resolution(pp);
	    else helpq=Majaxis(pp);
	    helpp+=helpq;
	    helpp/=2.;
	    DISTV1(helph,Pos(p),Pos(pp));
	    if (helph > helpp) pll=1;
	    pr=(bodyptr) Brothers(pq)[k];
	    for (iz=0;iz< Nbrothers(p);iz++) {
	      pp=(bodyptr) Brothers(p)[iz];
	      if (Majaxis(pr) == 0) helpp=Resolution(pr);
	      else helpp=Majaxis(pr);
	      if (Majaxis(pp) == 0) helpq=Resolution(pp);
	      else helpq=Majaxis(pp);
	      helpp+=helpq;
	      helpp/=2.;
	      DISTV1(helph,Pos(pr),Pos(pp));
	      if (helph > helpp) pll=1;
	    }
	    if (pll == 0) {
	      DISTV1(helph,Pos(p),Pos(Brothers(pq)[k]));
	      if (helph <= (Resolution(p)+Resolution(Brothers(pq)[k]))/2.) {
		Brothers(p)[Nbrothers(p)]=Brothers(pq)[k];
		Nbrothers(p)++;
	      }
	    }
	  }
	}
	flag[ii]=-1;
      }
      for (ii=0;ii<Nbrothers(p);ii++) {
	pq=(bodyptr) Brothers(p)[ii];
	if (Nbrothers(pq) > 0 && flag[ii] >= 0) {
	  numin1=Frequency(p);
	  for (k=0;k<Nbrothers(p);k++) numin1=MIN(numin1,Frequency(Brothers(p)[k]));
	  numax1=Frequency(p);
	  for (k=0;k<Nbrothers(p);k++) numax1=MAX(numax1,Frequency(Brothers(p)[k]));
	  numin2=Frequency(pq);
	  for (k=0;k<Nbrothers(pq);k++) numin2=MIN(numin2,Frequency(Brothers(pq)[k]));
	  numax2=Frequency(pq);
	  for (k=0;k<Nbrothers(pq);k++) numax2=MAX(numax2,Frequency(Brothers(pq)[k]));
	  if (numax1 < numin2*1.2 || numin1 > numax2*0.8 && fabs(Specindx(p)-Specindx(pq)) < 2.0) flag[ii]=1;
	}
      }
      for (ii=0;ii<Nbrothers(p);ii++) if (flag[ii] == 1) {
	pq=(bodyptr) Brothers(p)[ii];
	for (k=0;k<Nbrothers(pq);k++) {
	  nn=0;
	  strc=strcmp(Name(p),Name(Brothers(pq)[k]));
	  if (strc == 0) nn=1; 
	  for (ij=0;ij<Nbrothers(p);ij++) {
	    strc=strcmp(Name(Brothers(p)[ij]),Name(Brothers(pq)[k]));
	    if (strc == 0) nn=1;
	  }
	  if (nn == 0) {
	    pll=0;
	    pp=(bodyptr) Brothers(pq)[k];
	    if (Majaxis(p) == 0) helpp=Resolution(p);
	    else helpp=Majaxis(p);
	    if (Majaxis(pp) == 0) helpq=Resolution(pp);
	    else helpq=Majaxis(pp);
	    helpp+=helpq;
	    helpp/=2.;
	    DISTV1(helph,Pos(p),Pos(pp));
	    if (helph > helpp) pll=1;
	    pr=(bodyptr) Brothers(pq)[k];
	    for (iz=0;iz< Nbrothers(p);iz++) {
	      pp=(bodyptr) Brothers(p)[iz];
	      if (Majaxis(pr) == 0) helpp=Resolution(pr);
	      else helpp=Majaxis(pr);
	      if (Majaxis(pp) == 0) helpq=Resolution(pp);
	      else helpq=Majaxis(pp);
	      helpp+=helpq;
	      helpp/=2.;
	      DISTV1(helph,Pos(pr),Pos(pp));
	      if (helph > helpp) pll=1;
	    }
	    if (pll == 0) {
	      DISTV1(helph,Pos(p),Pos(Brothers(pq)[k]));
	      if (helph <= (Resolution(p)+Resolution(Brothers(pq)[k]))/2.) {
		Brothers(p)[Nbrothers(p)]=Brothers(pq)[k];
		Nbrothers(p)++;
	      }
	    }
	  }
	}
	flag[ii]=-1;
      }
      for (ii=0;ii<Nbrothers(p);ii++) {
	pq=(bodyptr) Brothers(p)[ii];
	if (Nbrothers(pq) > 0 && flag[ii] >= 0) {
	  jj=1;
	  freqhelp[0]=Frequency(p);
	  for (i=0;i<Nbrothers(p);i++) {
	    j=0;
	    for (kk=0;kk<jj;kk++) if (Frequency(Brothers(p)[i]) == freqhelp[kk]) j=1;
	    if (j == 0) {
	      freqhelp[jj]=Frequency(Brothers(p)[i]);
	      jj++;
	    }
	  }
	  ninit=jj;
	  jj=1;
	  freqhelp1[0]=Frequency(pq);
	  for (i=0;i<Nbrothers(pq);i++) {
	    j=0;
	    for (kk=0;kk<jj;kk++) if (Frequency(Brothers(pq)[i]) == freqhelp1[kk]) j=1;
	    if (j == 0) {
	      freqhelp1[jj]=Frequency(Brothers(pq)[i]);
	      jj++;
	    }
	  }    
	  nnfinal=jj;
	  if (ninit >= nnfinal && Specindx(p) != 0.0) {
	    totalspread++;
	    j=k=0;
	    for (k=0;k<Nbrothers(pq);k++) if (Frequency(Brothers(pq)[k]) == Frequency(pq)) j++;
	    pqq=pq;
	    jj=0;
	    for (i=0;i<Nbrothers(p);i++) {
	      strc=strcmp(Name(pq),Name(Brothers(p)[i]));
	      if (strc == 0) jj=1;
	    }
	    if (jj == 1) {
	      if (Frequency(p) != Frequency(pq)) {
		pll=0;
		pp=pq;
		if (Majaxis(p) == 0) helpp=Resolution(p);
		else helpp=Majaxis(p);
		if (Majaxis(pp) == 0) helpq=Resolution(pp);
		else helpq=Majaxis(pp);
		helpp+=helpq;
		helpp/=2.;
		DISTV1(helph,Pos(p),Pos(pp));
		if (helph > helpp) pll=1;
		pr=p;
		for (iz=0;iz< Nbrothers(pq);iz++) {
		  pp=(bodyptr) Brothers(pq)[iz];
		  if (Majaxis(pr) == 0) helpp=Resolution(pr);
		  else helpp=Majaxis(pr);
		  if (Majaxis(pp) == 0) helpq=Resolution(pp);
		  else helpq=Majaxis(pp);
		  helpp+=helpq;
		  helpp/=2.;
		  DISTV1(helph,Pos(pr),Pos(pp));
		  if (helph > helpp) pll=1;
		}
		if (pll == 0) {
		  DISTV1(helph,Pos(p),Pos(pq));
		  if (helph <= (Resolution(p)+Resolution(pq))/2.) {
		    Brothers(pq)[j]=(nodeptr) p;
		    j++;
		    k++;
		  }
		}
	      }
	      for (i=0;i<Nbrothers(p);i++) {
		if (Frequency(pq) != Frequency(Brothers(p)[i])) {
		  pll=0;
		  pr=(bodyptr) Brothers(p)[i];
		  pp=pq;
		  if (Majaxis(pr) == 0) helpp=Resolution(pr);
		  else helpp=Majaxis(pr);
		  if (Majaxis(pp) == 0) helpq=Resolution(pp);
		  else helpq=Majaxis(pp);
		  helpp+=helpq;
		  helpp/=2.;
		  DISTV1(helph,Pos(p),Pos(pp));
		  if (helph > helpp) pll=1;
		  pr=(bodyptr) Brothers(p)[i];
		  for (iz=0;iz< Nbrothers(pq);iz++) {
		    pp=(bodyptr) Brothers(pq)[iz];
		    if (Majaxis(pr) == 0) helpp=Resolution(pr);
		    else helpp=Majaxis(pr);
		    if (Majaxis(pp) == 0) helpq=Resolution(pp);
		    else helpq=Majaxis(pp);
		    helpp+=helpq;
		    helpp/=2.;
		    DISTV1(helph,Pos(pr),Pos(pp));
		    if (helph > helpp) pll=1;
		  }
		  if (pll == 0) {
		    DISTV1(helph,Pos(Brothers(p)[i]),Pos(pq));
		    if (helph <= (Resolution(Brothers(p)[i])+Resolution(pq))/2.) {
		      Brothers(pq)[j]=Brothers(p)[i];
		      j++;
		      k++;
		    }
		  }
		}
	      }
	      Nbrothers(pq)=j;
	      Specindx(pq)=Specindx(p);
	      Absz(pq)=Absz(p);
	      Resolved(pq)=Resolved(p);
	    }
	  } else if (ninit < nnfinal && Specindx(pq) != 0.0) {
	    totalspread++;
	    j=k=0;
	    for (k=0;k<Nbrothers(p);k++) if (Frequency(Brothers(p)[k]) == Frequency(p)) j++;
	    pqq=pq;
	    jj=0;
	    for (i=0;i<Nbrothers(pq);i++) {
	      strc=strcmp(Name(p),Name(Brothers(pq)[i]));
	      if (strc == 0) jj=1;
	    }
	    if (jj == 1) {
	      if (Frequency(pq) != Frequency(p)) {
		pll=0;
		pr=p;
		pp=pq;
		if (Majaxis(pr) == 0) helpp=Resolution(pr);
		else helpp=Majaxis(pr);
		if (Majaxis(pp) == 0) helpq=Resolution(pp);
		else helpq=Majaxis(pp);
		helpp+=helpq;
		helpp/=2.;
		DISTV1(helph,Pos(p),Pos(pp));
		if (helph > helpp) pll=1;
		pr=pq;
		for (iz=0;iz< Nbrothers(p);iz++) {
		  pp=(bodyptr) Brothers(p)[iz];
		  if (Majaxis(pr) == 0) helpp=Resolution(pr);
		  else helpp=Majaxis(pr);
		  if (Majaxis(pp) == 0) helpq=Resolution(pp);
		  else helpq=Majaxis(pp);
		  helpp+=helpq;
		  helpp/=2.;
		  DISTV1(helph,Pos(pr),Pos(pp));
		  if (helph > helpp) pll=1;
		}
		if (pll == 0) {
		  DISTV1(helph,Pos(p),Pos(pq));
		  if (helph <= (Resolution(p)+Resolution(pq))/2.) {
		    Brothers(p)[j]=(nodeptr) pq;
		    j++;
		    k++;
		  }
		}
	      }
	      for (i=0;i<Nbrothers(pq);i++) {
		if (Frequency(p) != Frequency(Brothers(pq)[i])) {
		  pll=0;
		  pr=(bodyptr) Brothers(pq)[i];
		  pp=p;
		  if (Majaxis(pr) == 0) helpp=Resolution(pr);
		  else helpp=Majaxis(pr);
		  if (Majaxis(pp) == 0) helpq=Resolution(pp);
		  else helpq=Majaxis(pp);
		  helpp+=helpq;
		  helpp/=2.;
		  DISTV1(helph,Pos(p),Pos(pp));
		  if (helph > helpp) pll=1;
		  pr=(bodyptr) Brothers(pq)[i];
		  for (iz=0;iz< Nbrothers(p);iz++) {
		    pp=(bodyptr) Brothers(p)[iz];
		    if (Majaxis(pr) == 0) helpp=Resolution(pr);
		    else helpp=Majaxis(pr);
		    if (Majaxis(pp) == 0) helpq=Resolution(pp);
		    else helpq=Majaxis(pp);
		    helpp+=helpq;
		    helpp/=2.;
		    DISTV1(helph,Pos(pr),Pos(pp));
		    if (helph > helpp) pll=1;
		  }
		  if (pll == 0) {
		    DISTV1(helph,Pos(p),Pos(Brothers(pq)[i]));
		    if (helph <= (Resolution(p)+Resolution(Brothers(pq)[i]))/2.) {
		      Brothers(p)[j]=Brothers(pq)[i];
		      j++;
		      k++;
		    }
		  }
		}
	      }
	      Nbrothers(p)=j;
	      Specindx(p)=Specindx(pq);
	      Absz(p)=Absz(pq);
	      Resolved(p)=Resolved(pq);
	    }
	  } else {
	    qwe=-1;
	    if (Specindx(pq) == 0) qwe=0;
	    else if (Specindx(p) == 0) qwe=1;
	    else {
	      if (Specindx(p) < 0.0 && Specindx(pq) < 0.0 && Specindx(p) > Specindx(pq)) qwe=0;
	      else if (Specindx(p) < 0.0 && Specindx(pq) < 0.0 && Specindx(p) < Specindx(pq)) qwe=1;
	      else if (Specindx(p)*Specindx(pq) < 0.0 && Specindx(p) < 0.0) qwe=0;
	      else if (Specindx(p)*Specindx(pq) < 0.0 && Specindx(pq) < 0.0) qwe=1;
	      else if (Specindx(p) > 0.0 && Specindx(pq) > 0.0 && Specindx(p) < Specindx(pq)) qwe=0;
	      else if (Specindx(p) > 0.0 && Specindx(pq) > 0.0 && Specindx(p) > Specindx(pq)) qwe=1;
	    }
	    if (qwe == -1) {
	      printf("problem Specindx(p)=%5.3e  Specindx(pq)=%5.3e \n",Specindx(p),Specindx(pq));
	      getchar();
	    }
	    if (qwe == 0) {
	      totalspread++;
	      j=k=0;
	      for (k=0;k<Nbrothers(pq);k++) if (Frequency(Brothers(pq)[k]) == Frequency(pq)) j++;
	      pqq=pq;
	      jj=0;
	      for (i=0;i<Nbrothers(p);i++) {
		strc=strcmp(Name(pq),Name(Brothers(p)[i]));
		if (strc == 0) jj=1;
	      }
	      if (jj == 1) {
		if (Frequency(p) != Frequency(pq)) {
		  pll=0;
		  pp=pq;
		  if (Majaxis(p) == 0) helpp=Resolution(p);
		  else helpp=Majaxis(p);
		  if (Majaxis(pp) == 0) helpq=Resolution(pp);
		  else helpq=Majaxis(pp);
		  helpp+=helpq;
		  helpp/=2.;
		  DISTV1(helph,Pos(p),Pos(pp));
		  if (helph > helpp) pll=1;
		  pr=p;
		  for (iz=0;iz< Nbrothers(pq);iz++) {
		    pp=(bodyptr) Brothers(pq)[iz];
		    if (Majaxis(pr) == 0) helpp=Resolution(pr);
		    else helpp=Majaxis(pr);
		    if (Majaxis(pp) == 0) helpq=Resolution(pp);
		    else helpq=Majaxis(pp);
		    helpp+=helpq;
		    helpp/=2.;
		    DISTV1(helph,Pos(pr),Pos(pp));
		    if (helph > helpp) pll=1;
		  }
		  if (pll == 0) {
		    DISTV1(helph,Pos(p),Pos(pq));
		    if (helph <= (Resolution(p)+Resolution(pq))/2.) {
		      Brothers(pq)[j]=(nodeptr) p;
		      j++;
		      k++;
		    }
		  }
		}
		for (i=0;i<Nbrothers(p);i++) {
		  if (Frequency(pq) != Frequency(Brothers(p)[i])) {
		    pll=0;
		    pr=(bodyptr) Brothers(p)[i];
		    pp=pq;
		    if (Majaxis(pr) == 0) helpp=Resolution(pr);
		    else helpp=Majaxis(pr);
		    if (Majaxis(pp) == 0) helpq=Resolution(pp);
		    else helpq=Majaxis(pp);
		    helpp+=helpq;
		    helpp/=2.;
		    DISTV1(helph,Pos(p),Pos(pp));
		    if (helph > helpp) pll=1;
		    pr=(bodyptr) Brothers(p)[i];
		    for (iz=0;iz< Nbrothers(pq);iz++) {
		      pp=(bodyptr) Brothers(pq)[iz];
		      if (Majaxis(pr) == 0) helpp=Resolution(pr);
		      else helpp=Majaxis(pr);
		      if (Majaxis(pp) == 0) helpq=Resolution(pp);
		      else helpq=Majaxis(pp);
		      helpp+=helpq;
		      helpp/=2.;
		      DISTV1(helph,Pos(pr),Pos(pp));
		      if (helph > helpp) pll=1;
		    }
		    if (pll == 0) {
		      DISTV1(helph,Pos(Brothers(p)[i]),Pos(pq));
		      if (helph <= (Resolution(Brothers(p)[i])+Resolution(pq))/2.) {
			Brothers(pq)[j]=Brothers(p)[i];
			j++;
			k++;
		      }
		    }
		  }
		}
		Nbrothers(pq)=j;
		Specindx(pq)=Specindx(p);
		Absz(pq)=Absz(p);
		Resolved(pq)=Resolved(p);
	      }
	    } else {
	      totalspread++;
	      j=k=0;
	      for (k=0;k<Nbrothers(p);k++) if (Frequency(Brothers(p)[k]) == Frequency(p)) j++;
	      pqq=pq;
	      jj=0;
	      for (i=0;i<Nbrothers(pq);i++) {
		strc=strcmp(Name(p),Name(Brothers(pq)[i]));
		if (strc == 0) jj=1;
	      }
	      if (jj == 1) {
		if (Frequency(pq) != Frequency(p)) {
		  pll=0;
		  pr=p;
		  pp=pq;
		  if (Majaxis(pr) == 0) helpp=Resolution(pr);
		  else helpp=Majaxis(pr);
		  if (Majaxis(pp) == 0) helpq=Resolution(pp);
		  else helpq=Majaxis(pp);
		  helpp+=helpq;
		  helpp/=2.;
		  DISTV1(helph,Pos(p),Pos(pp));
		  if (helph > helpp) pll=1;
		  pr=pq;
		  for (iz=0;iz< Nbrothers(p);iz++) {
		    pp=(bodyptr) Brothers(p)[iz];
		    if (Majaxis(pr) == 0) helpp=Resolution(pr);
		    else helpp=Majaxis(pr);
		    if (Majaxis(pp) == 0) helpq=Resolution(pp);
		    else helpq=Majaxis(pp);
		    helpp+=helpq;
		    helpp/=2.;
		    DISTV1(helph,Pos(pr),Pos(pp));
		    if (helph > helpp) pll=1;
		  }
		  if (pll == 0) {
		    DISTV1(helph,Pos(p),Pos(pq));
		    if (helph <= (Resolution(p)+Resolution(pq))/2.) {
		      Brothers(p)[j]=(nodeptr) pq;
		      j++;
		      k++;
		    }
		  }
		}
		for (i=0;i<Nbrothers(pq);i++) {
		  if (Frequency(p) != Frequency(Brothers(pq)[i])) {
		    pll=0;
		    pr=(bodyptr) Brothers(pq)[i];
		    pp=p;
		    if (Majaxis(pr) == 0) helpp=Resolution(pr);
		    else helpp=Majaxis(pr);
		    if (Majaxis(pp) == 0) helpq=Resolution(pp);
		    else helpq=Majaxis(pp);
		    helpp+=helpq;
		    helpp/=2.;
		    DISTV1(helph,Pos(p),Pos(pp));
		    if (helph > helpp) pll=1;
		    pr=(bodyptr) Brothers(pq)[i];
		    for (iz=0;iz< Nbrothers(p);iz++) {
		      pp=(bodyptr) Brothers(p)[iz];
		      if (Majaxis(pr) == 0) helpp=Resolution(pr);
		      else helpp=Majaxis(pr);
		      if (Majaxis(pp) == 0) helpq=Resolution(pp);
		      else helpq=Majaxis(pp);
		      helpp+=helpq;
		      helpp/=2.;
		      DISTV1(helph,Pos(pr),Pos(pp));
		      if (helph > helpp) pll=1;
		    }
		    if (pll == 0) {
		      DISTV1(helph,Pos(p),Pos(Brothers(pq)[i]));
		      if (helph <= (Resolution(p)+Resolution(Brothers(pq)[i]))/2.) {
			Brothers(p)[j]=Brothers(pq)[i];
			j++;
			k++;
		      }
		    }
		  }
		}
		Nbrothers(p)=j;
		Specindx(p)=Specindx(pq);
		Absz(p)=Absz(pq);
		Resolved(p)=Resolved(pq);
	      }
	    }    
	  }  
	}
      }
    }
  }
  
  for (p=bodytab;p<bodytab+nbody;p++) {  /* check for sources that are not brothers at the same frequency, but share the same spectrum (error bars) */
    if (Nbrothers(p) > 0) {
      for (ii=0;ii<Nbrothers(p);ii++) flag[ii]=0;
      for (ii=0;ii<Nbrothers(p);ii++) {
	pq=(bodyptr) Brothers(p)[ii];
	if (Nbrothers(pq) > 0) {
	  if (Specindx(p) > Specindx(pq)-0.3 && Specindx(p) < Specindx(pq)+0.3) flag[ii]=1;
	}
      }
      for (ii=0;ii<Nbrothers(p);ii++) if (flag[ii] == 1) {
	pq=(bodyptr) Brothers(p)[ii];
	nn=0;
	for (k=0;k<Nbrothers(pq);k++) {
	  strc=strcmp(Name(p),Name(Brothers(pq)[k]));
	  if (strc == 0) nn=1; 
	}
	if (nn == 0) { 
	  pll=0;
	  pr=p;
	  pp=pq;
	  if (Majaxis(pr) == 0) helpp=Resolution(pr);
	  else helpp=Majaxis(pr);
	  if (Majaxis(pp) == 0) helpq=Resolution(pp);
	  else helpq=Majaxis(pp);
	  helpp+=helpq;
	  helpp/=2.;
	  DISTV1(helph,Pos(p),Pos(pp));
	  if (helph > helpp) pll=1;
	  pr=p;
	  for (iz=0;iz< Nbrothers(pq);iz++) {
	    pp=(bodyptr) Brothers(pq)[iz];
	    if (Majaxis(pr) == 0) helpp=Resolution(pr);
	    else helpp=Majaxis(pr);
	    if (Majaxis(pp) == 0) helpq=Resolution(pp);
	    else helpq=Majaxis(pp);
	    helpp+=helpq;
	    helpp/=2.;
	    DISTV1(helph,Pos(pr),Pos(pp));
	    if (helph > helpp) pll=1;
	  }
	  if (pll == 0) {
	    DISTV1(helph,Pos(p),Pos(pq));
	    if (helph <= (Resolution(p)+Resolution(pq))/2.) {
	      Brothers(pq)[Nbrothers(pq)]=(nodeptr) p;
	      Nbrothers(pq)++;
	    }
	  }
	} 
	flag[ii]=-1;
      }
      for (ii=0;ii<Nbrothers(p);ii++) {
	pq=(bodyptr) Brothers(p)[ii];
	if (Nbrothers(pq) > 0 && flag[ii] >= 0) {
	  numin1=Frequency(p);
	  for (k=0;k<Nbrothers(p);k++) numin1=MIN(numin1,Frequency(Brothers(p)[k]));
	  numax1=Frequency(p);
	  for (k=0;k<Nbrothers(p);k++) numax1=MAX(numax1,Frequency(Brothers(p)[k]));
	  numin2=Frequency(pq);
	  for (k=0;k<Nbrothers(pq);k++) numin2=MIN(numin2,Frequency(Brothers(pq)[k]));
	  numax2=Frequency(pq);
	  for (k=0;k<Nbrothers(pq);k++) numax2=MAX(numax2,Frequency(Brothers(pq)[k]));
	  if (numax1 < numin2*1.2 || numin1 > numax2*0.8 && fabs(Specindx(p)-Specindx(pq)) < 2.0) flag[ii]=1;
	}
      }
      for (ii=0;ii<Nbrothers(p);ii++) if (flag[ii] == 1) {
	pq=(bodyptr) Brothers(p)[ii];
	nn=0;
	for (k=0;k<Nbrothers(pq);k++) {
	  strc=strcmp(Name(p),Name(Brothers(pq)[k]));
	  if (strc == 0) nn=1; 
	}
	if (nn == 0) { 
	  pll=0;
	  pr=p;
	  pp=pq;
	  if (Majaxis(pr) == 0) helpp=Resolution(pr);
	  else helpp=Majaxis(pr);
	  if (Majaxis(pp) == 0) helpq=Resolution(pp);
	  else helpq=Majaxis(pp);
	  helpp+=helpq;
	  helpp/=2.;
	  DISTV1(helph,Pos(p),Pos(pp));
	  if (helph > helpp) pll=1;
	  pr=p;
	  for (iz=0;iz< Nbrothers(pq);iz++) {
	    pp=(bodyptr) Brothers(pq)[iz];
	    if (Majaxis(pr) == 0) helpp=Resolution(pr);
	    else helpp=Majaxis(pr);
	    if (Majaxis(pp) == 0) helpq=Resolution(pp);
	    else helpq=Majaxis(pp);
	    helpp+=helpq;
	    helpp/=2.;
	    DISTV1(helph,Pos(pr),Pos(pp));
	    if (helph > helpp) pll=1;
	  }
	  if (pll == 0) {
	    DISTV1(helph,Pos(p),Pos(pq));
	    if (helph <= (Resolution(p)+Resolution(pq))/2.) {
	      Brothers(pq)[Nbrothers(pq)]=(nodeptr) p;
	      Nbrothers(pq)++;
	    }
	  }
	} 
	flag[ii]=-1;
      }
      for (ii=0;ii<Nbrothers(p);ii++) {
	pq=(bodyptr) Brothers(p)[ii];
	if (Nbrothers(pq) > 0 && flag[ii] >= 0) {
	  jj=1;
	  freqhelp[0]=Frequency(p);
	  for (i=0;i<Nbrothers(p);i++) {
	    j=0;
	    for (kk=0;kk<jj;kk++) if (Frequency(Brothers(p)[i]) == freqhelp[kk]) j=1;
	    if (j == 0) {
	      freqhelp[jj]=Frequency(Brothers(p)[i]);
	      jj++;
	    }
	  }
	  ninit=jj;
	  jj=1;
	  freqhelp1[0]=Frequency(pq);
	  for (i=0;i<Nbrothers(pq);i++) {
	    j=0;
	    for (kk=0;kk<jj;kk++) if (Frequency(Brothers(pq)[i]) == freqhelp1[kk]) j=1;
	    if (j == 0) {
	      freqhelp1[jj]=Frequency(Brothers(pq)[i]);
	      jj++;
	    }
	  }    
	  nnfinal=jj;
	  if (ninit >= nnfinal && Specindx(p) != 0.0) {
	    totalspread++;
	    j=k=0;
	    for (k=0;k<Nbrothers(pq);k++) if (Frequency(Brothers(pq)[k]) == Frequency(pq)) j++;
	    pqq=pq;
	    jj=0;
	    for (i=0;i<Nbrothers(p);i++) {
	      strc=strcmp(Name(pq),Name(Brothers(p)[i]));
	      if (strc == 0) jj=1;
	    }
	    if (jj == 1) {
	      if (Frequency(p) != Frequency(pq)) {
		pll=0;
		pr=p;
		pp=pq;
		if (Majaxis(pr) == 0) helpp=Resolution(pr);
		else helpp=Majaxis(pr);
		if (Majaxis(pp) == 0) helpq=Resolution(pp);
		else helpq=Majaxis(pp);
		helpp+=helpq;
		helpp/=2.;
		DISTV1(helph,Pos(p),Pos(pp));
		if (helph > helpp) pll=1;
		pr=p;
		for (iz=0;iz< Nbrothers(pq);iz++) {
		  pp=(bodyptr) Brothers(pq)[iz];
		  if (Majaxis(pr) == 0) helpp=Resolution(pr);
		  else helpp=Majaxis(pr);
		  if (Majaxis(pp) == 0) helpq=Resolution(pp);
		  else helpq=Majaxis(pp);
		  helpp+=helpq;
		  helpp/=2.;
		  DISTV1(helph,Pos(pr),Pos(pp));
		  if (helph > helpp) pll=1;
		}
		if (pll == 0) {
		  DISTV1(helph,Pos(p),Pos(pq));
		  if (helph <= (Resolution(p)+Resolution(pq))/2.) {
		    Brothers(pq)[j]=(nodeptr) p;
		    j++;
		    k++;
		  }
		}
	      }
	      for (i=0;i<Nbrothers(p);i++) {
		if (Frequency(pq) != Frequency(Brothers(p)[i])) {
		  pll=0;
		  pr=(bodyptr) Brothers(p)[i];
		  pp=pq;
		  if (Majaxis(pr) == 0) helpp=Resolution(pr);
		  else helpp=Majaxis(pr);
		  if (Majaxis(pp) == 0) helpq=Resolution(pp);
		  else helpq=Majaxis(pp);
		  helpp+=helpq;
		  helpp/=2.;
		  DISTV1(helph,Pos(p),Pos(pp));
		  if (helph > helpp) pll=1;
		  pr=(bodyptr) Brothers(p)[i];
		  for (iz=0;iz< Nbrothers(pq);iz++) {
		    pp=(bodyptr) Brothers(pq)[iz];
		    if (Majaxis(pr) == 0) helpp=Resolution(pr);
		    else helpp=Majaxis(pr);
		    if (Majaxis(pp) == 0) helpq=Resolution(pp);
		    else helpq=Majaxis(pp);
		    helpp+=helpq;
		    helpp/=2.;
		    DISTV1(helph,Pos(pr),Pos(pp));
		    if (helph > helpp) pll=1;
		  }
		  if (pll == 0) {
		    DISTV1(helph,Pos(pq),Pos(Brothers(p)[i]));
		    if (helph <= (Resolution(pq)+Resolution(Brothers(p)[i]))/2.) {
		      Brothers(pq)[j]=Brothers(p)[i];
		      j++;
		      k++;
		    }
		  }
		}
	      }
	      Nbrothers(pq)=j;
	      Specindx(pq)=Specindx(p);
	      Absz(pq)=Absz(p);
	      Resolved(pq)=Resolved(p);
	    }
	  } else if (ninit < nnfinal && Specindx(pq) != 0.0) {
	    totalspread++;
	    j=k=0;
	    for (k=0;k<Nbrothers(p);k++) if (Frequency(Brothers(p)[k]) == Frequency(p)) j++;
	    pqq=pq;
	    jj=0;
	    for (i=0;i<Nbrothers(pq);i++) {
	      strc=strcmp(Name(p),Name(Brothers(pq)[i]));
	      if (strc == 0) jj=1;
	    }
	    if (jj == 1) {
	      if (Frequency(pq) != Frequency(p)) {
		pll=0;
		pr=pq;
		pp=p;
		if (Majaxis(pr) == 0) helpp=Resolution(pr);
		else helpp=Majaxis(pr);
		if (Majaxis(pp) == 0) helpq=Resolution(pp);
		else helpq=Majaxis(pp);
		helpp+=helpq;
		helpp/=2.;
		DISTV1(helph,Pos(p),Pos(pp));
		if (helph > helpp) pll=1;
		pr=pq;
		for (iz=0;iz< Nbrothers(p);iz++) {
		  pp=(bodyptr) Brothers(p)[iz];
		  if (Majaxis(pr) == 0) helpp=Resolution(pr);
		  else helpp=Majaxis(pr);
		  if (Majaxis(pp) == 0) helpq=Resolution(pp);
		  else helpq=Majaxis(pp);
		  helpp+=helpq;
		  helpp/=2.;
		  DISTV1(helph,Pos(pr),Pos(pp));
		  if (helph > helpp) pll=1;
		}
		if (pll == 0) {
		  DISTV1(helph,Pos(p),Pos(pq));
		  if (helph <= (Resolution(p)+Resolution(pq))/2.) {
		    Brothers(p)[j]=(nodeptr) pq;
		    j++;
		    k++;
		  }
		}
	      }
	      for (i=0;i<Nbrothers(pq);i++) {
		if (Frequency(p) != Frequency(Brothers(pq)[i])) {
		  pll=0;
		  pr=(bodyptr) Brothers(pq)[i];
		  pp=p;
		  if (Majaxis(pr) == 0) helpp=Resolution(pr);
		  else helpp=Majaxis(pr);
		  if (Majaxis(pp) == 0) helpq=Resolution(pp);
		  else helpq=Majaxis(pp);
		  helpp+=helpq;
		  helpp/=2.;
		  DISTV1(helph,Pos(p),Pos(pp));
		  if (helph > helpp) pll=1;
		  pr=(bodyptr) Brothers(pq)[i];
		  for (iz=0;iz< Nbrothers(p);iz++) {
		    pp=(bodyptr) Brothers(p)[iz];
		    if (Majaxis(pr) == 0) helpp=Resolution(pr);
		    else helpp=Majaxis(pr);
		    if (Majaxis(pp) == 0) helpq=Resolution(pp);
		    else helpq=Majaxis(pp);
		    helpp+=helpq;
		    helpp/=2.;
		    DISTV1(helph,Pos(pr),Pos(pp));
		    if (helph > helpp) pll=1;
		  }
		  if (pll == 0) {
		    DISTV1(helph,Pos(p),Pos(Brothers(pq)[i]));
		    if (helph <= (Resolution(p)+Resolution(Brothers(pq)[i]))/2.) {
		      Brothers(p)[j]=Brothers(pq)[i];
		      j++;
		      k++;
		    }
		  }
		}
	      }
	      Nbrothers(p)=j;
	      Specindx(p)=Specindx(pq);
	      Absz(p)=Absz(pq);
	      Resolved(p)=Resolved(pq);
	    }
	  } else {
	    qwe=-1;
	    if (Specindx(pq) == 0) qwe=0;
	    else if (Specindx(p) == 0) qwe=1;
	    else {
	      if (Specindx(p) < 0.0 && Specindx(pq) < 0.0 && Specindx(p) > Specindx(pq)) qwe=0;
	      else if (Specindx(p) < 0.0 && Specindx(pq) < 0.0 && Specindx(p) < Specindx(pq)) qwe=1;
	      else if (Specindx(p)*Specindx(pq) < 0.0 && Specindx(p) < 0.0) qwe=0;
	      else if (Specindx(p)*Specindx(pq) < 0.0 && Specindx(pq) < 0.0) qwe=1;
	      else if (Specindx(p) > 0.0 && Specindx(pq) > 0.0 && Specindx(p) < Specindx(pq)) qwe=0;
	      else if (Specindx(p) > 0.0 && Specindx(pq) > 0.0 && Specindx(p) > Specindx(pq)) qwe=1;
	    }
	    if (qwe == -1) {
	      printf("problem1 Specindx(p)=%5.3e  Specindx(pq)=%5.3e \n",Specindx(p),Specindx(pq));
	      getchar();
	    }
	    if (qwe == 0) {
	      totalspread++;
	      j=k=0;
	      for (k=0;k<Nbrothers(pq);k++) if (Frequency(Brothers(pq)[k]) == Frequency(pq)) j++;
	      pqq=pq;
	      jj=0;
	      for (i=0;i<Nbrothers(p);i++) {
		strc=strcmp(Name(pq),Name(Brothers(p)[i]));
		if (strc == 0) jj=1;
	      }
	      if (jj == 1) {
		if (Frequency(p) != Frequency(pq)) {
		  pll=0;
		  pr=p;
		  pp=pq;
		  if (Majaxis(pr) == 0) helpp=Resolution(pr);
		  else helpp=Majaxis(pr);
		  if (Majaxis(pp) == 0) helpq=Resolution(pp);
		  else helpq=Majaxis(pp);
		  helpp+=helpq;
		  helpp/=2.;
		  DISTV1(helph,Pos(p),Pos(pp));
		  if (helph > helpp) pll=1;
		  pr=p;
		  for (iz=0;iz< Nbrothers(pq);iz++) {
		    pp=(bodyptr) Brothers(pq)[iz];
		    if (Majaxis(pr) == 0) helpp=Resolution(pr);
		    else helpp=Majaxis(pr);
		    if (Majaxis(pp) == 0) helpq=Resolution(pp);
		    else helpq=Majaxis(pp);
		    helpp+=helpq;
		    helpp/=2.;
		    DISTV1(helph,Pos(pr),Pos(pp));
		    if (helph > helpp) pll=1;
		  }
		  if (pll == 0) {
		    DISTV1(helph,Pos(p),Pos(pq));
		    if (helph <= (Resolution(p)+Resolution(pq))/2.) {
		      Brothers(pq)[j]=(nodeptr) p;
		      j++;
		      k++;
		    }
		  }
		}
		for (i=0;i<Nbrothers(p);i++) {
		  if (Frequency(pq) != Frequency(Brothers(p)[i])) {
		    pll=0;
		    pr=(bodyptr) Brothers(p)[i];
		    pp=pq;
		    if (Majaxis(pr) == 0) helpp=Resolution(pr);
		    else helpp=Majaxis(pr);
		    if (Majaxis(pp) == 0) helpq=Resolution(pp);
		    else helpq=Majaxis(pp);
		    helpp+=helpq;
		    helpp/=2.;
		    DISTV1(helph,Pos(p),Pos(pp));
		    if (helph > helpp) pll=1;
		    pr=(bodyptr) Brothers(p)[i];
		    for (iz=0;iz< Nbrothers(pq);iz++) {
		      pp=(bodyptr) Brothers(pq)[iz];
		      if (Majaxis(pr) == 0) helpp=Resolution(pr);
		      else helpp=Majaxis(pr);
		      if (Majaxis(pp) == 0) helpq=Resolution(pp);
		      else helpq=Majaxis(pp);
		      helpp+=helpq;
		      helpp/=2.;
		      DISTV1(helph,Pos(pr),Pos(pp));
		      if (helph > helpp) pll=1;
		    }
		    if (pll == 0) {
		      DISTV1(helph,Pos(pq),Pos(Brothers(p)[i]));
		      if (helph <= (Resolution(pq)+Resolution(Brothers(p)[i]))/2.) {
			Brothers(pq)[j]=Brothers(p)[i];
			j++;
			k++;
		      }
		    }
		  }
		}
		Nbrothers(pq)=j;
		Specindx(pq)=Specindx(p);
		Absz(pq)=Absz(p);
		Resolved(pq)=Resolved(p);
	      }
	    } else {
	      totalspread++;
	      j=k=0;
	      for (k=0;k<Nbrothers(p);k++) if (Frequency(Brothers(p)[k]) == Frequency(p)) j++;
	      pqq=pq;
	      jj=0;
	      for (i=0;i<Nbrothers(pq);i++) {
		strc=strcmp(Name(p),Name(Brothers(pq)[i]));
		if (strc == 0) jj=1;
	      }
	      if (jj == 1) {
		if (Frequency(pq) != Frequency(p)) {
		  pll=0;
		  pr=pq;
		  pp=p;
		  if (Majaxis(pr) == 0) helpp=Resolution(pr);
		  else helpp=Majaxis(pr);
		  if (Majaxis(pp) == 0) helpq=Resolution(pp);
		  else helpq=Majaxis(pp);
		  helpp+=helpq;
		  helpp/=2.;
		  DISTV1(helph,Pos(p),Pos(pp));
		  if (helph > helpp) pll=1;
		  pr=pq;
		  for (iz=0;iz< Nbrothers(p);iz++) {
		    pp=(bodyptr) Brothers(p)[iz];
		    if (Majaxis(pr) == 0) helpp=Resolution(pr);
		    else helpp=Majaxis(pr);
		    if (Majaxis(pp) == 0) helpq=Resolution(pp);
		    else helpq=Majaxis(pp);
		    helpp+=helpq;
		    helpp/=2.;
		    DISTV1(helph,Pos(pr),Pos(pp));
		    if (helph > helpp) pll=1;
		  }
		  if (pll == 0) {
		    DISTV1(helph,Pos(p),Pos(pq));
		    if (helph <= (Resolution(p)+Resolution(pq))/2.) {
		      Brothers(p)[j]=(nodeptr) pq;
		      j++;
		      k++;
		    }
		  }
		}
		for (i=0;i<Nbrothers(pq);i++) {
		  if (Frequency(p) != Frequency(Brothers(pq)[i])) {
		    pll=0;
		    pr=(bodyptr) Brothers(pq)[i];
		    pp=p;
		    if (Majaxis(pr) == 0) helpp=Resolution(pr);
		    else helpp=Majaxis(pr);
		    if (Majaxis(pp) == 0) helpq=Resolution(pp);
		    else helpq=Majaxis(pp);
		    helpp+=helpq;
		    helpp/=2.;
		    DISTV1(helph,Pos(p),Pos(pp));
		    if (helph > helpp) pll=1;
		    pr=(bodyptr) Brothers(pq)[i];
		    for (iz=0;iz< Nbrothers(p);iz++) {
		      pp=(bodyptr) Brothers(p)[iz];
		      if (Majaxis(pr) == 0) helpp=Resolution(pr);
		      else helpp=Majaxis(pr);
		      if (Majaxis(pp) == 0) helpq=Resolution(pp);
		      else helpq=Majaxis(pp);
		      helpp+=helpq;
		      helpp/=2.;
		      DISTV1(helph,Pos(pr),Pos(pp));
		      if (helph > helpp) pll=1;
		    }
		    if (pll == 0) {
		      DISTV1(helph,Pos(p),Pos(Brothers(pq)[i]));
		      if (helph <= (Resolution(p)+Resolution(Brothers(pq)[i]))/2.) {
			Brothers(p)[j]=Brothers(pq)[i];
			j++;
			k++;
		      }
		    }
		  }
		}
		Nbrothers(p)=j;
		Specindx(p)=Specindx(pq);
		Absz(p)=Absz(pq);
		Resolved(p)=Resolved(pq);
	      }
	    }    
	  }  
	}
      }
    }
  }

  for (p=bodytab;p<bodytab+nbody;p++) {	  
    if (Nbrothers(p) > 0 && Specindx(p) != 0.0) {
      helpp=Absz(p)+Specindx(p)*log10(Frequency(p));
      helpp=pow(10.,helpp);
      help=fabs(Flux(p)-helpp);
      if (help > 5.*Fluxerror(p)) {
	jj=0;
	while (Frequency(p) == Frequency(Brothers(p)[jj]) && jj < Nbrothers(p)-1) jj++;
	Nbrothers(p)=jj;
	Specindx(p)=0.0;
	Absz(p)=0.0;
      }
    }
  }  

  for (p=bodytab;p<bodytab+nbody;p++) {
    nn=0;
    if (Nbrothers(p) > 0) {
      for (i=0;i<Nbrothers(p);i++) {
	init[i]=(bodyptr) Brothers(p)[i];
	flag[i]=1;
      }
      for (i=0;i<Nbrothers(p);i++) {
	pq=(bodyptr) Brothers(p)[i];	  
	helpp=Absz(p)+Specindx(p)*log10(Frequency(pq));
	helpp=pow(10.,helpp);
	help=fabs(Flux(pq)-helpp);
	if (help > 5.*Fluxerror(pq)) {
	  flag[i]=0;
	  nn=1;
	}
      }
      if (nn == 1) {
	j=0;
	for (i=0;i<Nbrothers(p);i++) {
	  if (flag[i] == 1) {
	    Brothers(p)[j]=(nodeptr) init[i];
	    j++;
	  }
	}
	Nbrothers(p)=j;
      }
    }
  }	    

  for (p=bodytab;p<bodytab+nbody;p++) twofreq(p); /* security: remove source with only two different frequencies */
  
  for (p=bodytab;p<bodytab+nbody;p++) {
    if (Nbrothers(p) > 1) {
      if (Majaxis(p) == 0) helph=Resolution(p);
      else helph=Majaxis(p);
      npoints=1;
      for (i=0;i<Nbrothers(p);i++) {
	if (Frequency(p) != Frequency(Brothers(p)[i])) {
	  npoints++;
	  if (Majaxis(Brothers(p)[i]) == 0) helph=Resolution(Brothers(p)[i]);
	  else helph=Majaxis(Brothers(p)[i]);
	  helpp=MAX(helpp,helph);
	  helpq=MIN(helpq,helph);
	}
      }
      if (helpp/helpq > 10. && Specindx(p) < -2.0 && npoints == 3) {
	jj=0;
	while (Frequency(p) == Frequency(Brothers(p)[jj]) && jj < Nbrothers(p)-1) jj++;
	Nbrothers(p)=jj;
	Specindx(p)=0.0;
	Absz(p)=0.0;
      }
    }
  }

  for (p=bodytab;p<bodytab+nbody;p++) { /* remove double sources with the same name */
    if (Nbrothers(p) > 0) {
      for (i=0;i<51;i++) flag[i]=1;
      init[0]=p;
      for (k=0;k<Nbrothers(p);k++) init[k+1]=(bodyptr) Brothers(p)[k];
      for (k=0;k<Nbrothers(p);k++) {
	for (i=k+1;i<=Nbrothers(p);i++) {
	  strc=strcmp(Name(init[k]),Name(init[i]));
	  if (strc == 0 && flag[i] == 1) flag[i]=0;
	}
      }
      j=0;
      for (k=1;k<=Nbrothers(p);k++) if (flag[k] == 1) {
	Brothers(p)[j]=(nodeptr) init[k];
	j++;
      }
      Nbrothers(p)=j;
    }
  }

  doing=0;
  while (doing == 0) {
    doing=1;
    for (p=bodytab;p<bodytab+nbody;p++) {
      if (Nbrothers(p) > 0) {
	for (i=0;i<Nbrothers(p);i++) {
	  pq=(bodyptr) Brothers(p)[i];	  
	  if (Nbrothers(pq) > 0) {
	    nn=0;
	    for (k=0;k<Nbrothers(pq);k++) {
	      strc=strcmp(Name(Brothers(pq)[k]),Name(p));
	      if (strc == 0) nn=1;
	      for (ii=0;ii<Nbrothers(p);ii++) {
		strc=strcmp(Name(Brothers(pq)[k]),Name(Brothers(p)[ii]));
		if (strc == 0) nn=1;
	      }
	    }
	    if (nn == 0) {
	      jj=1;
	      freqhelp[0]=Frequency(p);
	      for (i=0;i<Nbrothers(p);i++) {
	      j=0;
	      for (kk=0;kk<jj;kk++) if (Frequency(Brothers(p)[i]) == freqhelp[kk]) j=1;
	      if (j == 0) {
		freqhelp[jj]=Frequency(Brothers(p)[i]);
		jj++;
	      }
	      }
	      ninit=jj;
	      jj=1;
	      freqhelp1[0]=Frequency(pq);
	      for (i=0;i<Nbrothers(pq);i++) {
		j=0;
		for (kk=0;kk<jj;kk++) if (Frequency(Brothers(pq)[i]) == freqhelp1[kk]) j=1;
		if (j == 0) {
		  freqhelp1[jj]=Frequency(Brothers(pq)[i]);
		  jj++;
		}
	      }    
	      nnfinal=jj;
	      qwe=-1;
	      if (Specindx(p) == 0.0) qwe=1;
	      else if (Specindx(pq) == 0.0) qwe=0;
	      else {
		if (ninit > nnfinal) qwe=0;
		else if (ninit < nnfinal) qwe=1;
		else {
		  if (Nbrothers(p) > Nbrothers(pq)) qwe=0;
		  else if (Nbrothers(p) < Nbrothers(pq)) qwe=1;
		  else {
		    if (fabs(Specindx(p)+0.9) <= fabs(Specindx(pq)+0.9)) qwe=0;
		    else qwe=1;
		  }
		}
	      }
	      if (qwe == -1) {
		printf("problem in spreadspec...\n");
		getchar();
	      }
	      if (qwe == 0) {
		for (ii=0;ii<Nbrothers(p);ii++) init[ii]=(bodyptr) Brothers(p)[ii];
		jj=0;
		for (ii=0;ii<Nbrothers(p);ii++) {
		  if (ii != i) {
		    Brothers(p)[jj]=(nodeptr) init[ii];
		    jj++;
		  }
		}
		Nbrothers(p)=jj;
		jj=0;
		while (Frequency(pq) == Frequency(Brothers(pq)[jj]) && jj < Nbrothers(pq)-1) jj++;
		Nbrothers(pq)=jj;
		Specindx(pq)=0.0;
		Absz(p)=0.0;
	      } else {
		for (ii=0;ii<Nbrothers(pq);ii++) init[ii]=(bodyptr) Brothers(pq)[ii];
		jj=0;
		for (ii=0;ii<Nbrothers(pq);ii++) {
		  strc=strcmp(Name(p),Name(Brothers(pq)[ii]));
		  if (strc != 0) Brothers(pq)[jj]=(nodeptr) init[ii];
		}
		Nbrothers(pq)=jj;
		jj=0;
		while (Frequency(p) == Frequency(Brothers(p)[jj]) && jj < Nbrothers(p)-1) jj++;
		Nbrothers(p)=jj;
		Specindx(p)=0.0;
		Absz(p)=0.0;
	      }
	    } 
	  }
	}
      }
    }
    
  
    for (p=bodytab;p<bodytab+nbody;p++) twofreq(p); /* security: remove source with only two different frequencies */
    
    for (p=bodytab;p<bodytab+nbody;p++) {
      if (Nbrothers(p) > 0) {
	for (i=0;i<Nbrothers(p);i++) {
	  pq=(bodyptr) Brothers(p)[i];
	  if (Nbrothers(pq) > 0) {
	    nn=0;
	    for (k=0;k<Nbrothers(pq);k++) {
	      strc=strcmp(Name(Brothers(pq)[k]),Name(p));
	      if (strc == 0) nn=1;
	      for (ii=0;ii<Nbrothers(p);ii++) {
		strc=strcmp(Name(Brothers(pq)[k]),Name(Brothers(p)[ii]));
		if (strc == 0) nn=1;
	      }
	    }
	    if (nn == 0) {
	      doing=0;
	    } 
	  }
	}
      }
    }
  }

  /* additional check - last cross check*/

  for (p=bodytab;p<bodytab+nbody;p++) {
    for (i=0;i<Nbrothers(p);i++) init[i]=(bodyptr) Brothers(p)[i];
    for (i=0;i<Nbrothers(p);i++) {
      if (Frequency(Brothers(p)[i]) != Frequency(p) && Specindx(Brothers(p)[i]) == 0.0) flag[i]=0;
      else flag[i]=1;
    }
    j=0;
    for (i=0;i<Nbrothers(p);i++) {
      if (flag[i] == 1) {
	Brothers(p)[j]=(nodeptr) init[i];
	j++;
      }
    }
    Nbrothers(p)=j;
    twofreq(p);
  }

  /*goto vorbei;*/

  printf("begin maketree with %i particles \n",nbody);
  excor=cmore=0;
  maketr=1;
  maketree(bodytab,nbody);       /* built up of the tree */
  maketr=0;

  printf("maketree done...\n");

  printf("last cross check \n");
  for (p=bodytab;p<bodytab+nbody;p++) {
    if (Specindx(p) != 0.0) {
      resold=Resolution(p);
      Resolution(p)=1200.;
      hackgrav1(p, Mass(p) > 0.0);    /* proximity search */	  
      Resolution(p)=resold;
      for (in=0;in<ccount;in++) {
	pq=(bodyptr) store[in];
	strc=strcmp(Name(p),Name(pq));
	if (strc != 0) {
	  if (Specindx(pq) != 0.0) {
	    strc1=0;
	    for (i=0;i<Nbrothers(p);i++) {
	      strc=strcmp(Name(Brothers(p)[i]),Name(pq));
	      if (strc == 0) strc1++;
	    }
	    for (j=0;j<Nbrothers(pq);j++) {
	      strc=strcmp(Name(p),Name(Brothers(pq)[j]));
	      if (strc == 0) strc1++;
	    }
	    for (i=0;i<Nbrothers(p);i++) {
	      for (j=0;j<Nbrothers(pq);j++) {
		strc=strcmp(Name(Brothers(p)[i]),Name(Brothers(pq)[j]));
		if (strc == 0) strc1++;
	      }
	    }
	    if (strc1 > 0) {
	      if (strc1 < Nbrothers(p)+1 || strc1 < Nbrothers(pq)+1) {
		if (fabs(Specindx(p)-Specindx(pq)) < 0.3) {		  
		  nn1=0; 
		  for (i=0;i<Nbrothers(pq);i++) {
		    strc=strcmp(Name(p),Name(Brothers(pq)[i]));
		    if (strc == 0) nn1++;
		  }
		  if (nn1 == 0) {
		    distance(p,pq);
		    if (diist == 0) {
		      Brothers(pq)[Nbrothers(pq)]=(nodeptr) p;
		      Nbrothers(pq)++;
		    }
		  }
		  nn1=0;
		  for (j=0;j<Nbrothers(p);j++) {
		    nn1=0;
		    strc=strcmp(Name(Brothers(p)[j]),Name(pq));
		    if (strc == 0) nn1++;
		    for (i=0;i<Nbrothers(pq);i++) {
		      strc=strcmp(Name(Brothers(p)[j]),Name(Brothers(pq)[i]));
		      if (strc == 0) nn1++;
		    }
		    if (nn1 == 0) {
		      distance((bodyptr) Brothers(p)[j],pq);
		      if (diist == 0) {
			Brothers(pq)[Nbrothers(pq)]=Brothers(p)[j];
			Nbrothers(pq)++;
		      }
		    }
		  }
		  nn1=0; 
		  for (i=0;i<Nbrothers(p);i++) {
		    strc=strcmp(Name(pq),Name(Brothers(p)[i]));
		    if (strc == 0) nn1++;
		  }
		  if (nn1 == 0) {
		    distance(pq, p);
		    if (diist == 0) {
		      Brothers(p)[Nbrothers(p)]=(nodeptr) pq;
		      Nbrothers(p)++;
		    }
		  }
		  nn1=0;
		  for (j=0;j<Nbrothers(pq);j++) {
		    nn1=0;
		    strc=strcmp(Name(Brothers(pq)[j]),Name(p));
		    if (strc == 0) nn1++;
		    for (i=0;i<Nbrothers(p);i++) {
		      strc=strcmp(Name(Brothers(pq)[j]),Name(Brothers(p)[i]));
		      if (strc == 0) nn1++;
		    }
		    if (nn1 == 0) {
		      distance((bodyptr) Brothers(pq)[j], p);
		      if (diist == 0) {
			Brothers(p)[Nbrothers(p)]=Brothers(pq)[j];
			Nbrothers(p)++;
		      }
		    }
		  }
		} else if (Nbrothers(p) < Nbrothers(pq)) {
		  for (i=0;i<Nbrothers(pq);i++) init[i]=(bodyptr) Brothers(pq)[i];
		  for (i=0;i<Nbrothers(pq);i++) {
		    strc=strcmp(Name(p),Name(Brothers(pq)[i]));
		    if (strc == 0) flag[i]=0;
		    else flag[i]=1;
		  }
		  j=0;
		  for (i=0;i<Nbrothers(pq);i++) {
		    if (flag[i] == 1) {
		      Brothers(pq)[j]=(nodeptr) init[i];
		      j++;
		    }
		  }
		  Nbrothers(pq)=j;
		  twofreq(pq);

		  for (i=0;i<Nbrothers(p);i++) init[i]=(bodyptr) Brothers(p)[i];
		  nn1=0;
		  strc=strcmp(Name(p),Name(pq));
		  if (strc == 0) nn1=1;
		  for (i=0;i<Nbrothers(pq);i++) {
		    strc=strcmp(Name(p),Name(Brothers(pq)[i]));
		    if (strc == 0) nn1=1;
		  }
		  if (nn1 == 1) {
		    for (i=0;i<Nbrothers(p);i++) if (Frequency(p) == Frequency(Brothers(p)[i])) flag[i]=1;
		    else flag[i]=0;
		    j=0;
		    for (i=0;i<Nbrothers(p);i++) {
		      if (flag[i] == 1) {
			Brothers(p)[j]=(nodeptr) init[i];
			j++;
		      }
		    }
		    Nbrothers(p)=j;
		    Specindx(p)=0.0;
		    Absz(p)=0.0;
		  } else {
		    for (i=0;i<Nbrothers(p);i++) init[i]=(bodyptr) Brothers(p)[i];
		    for (j=0;j<Nbrothers(p);j++) {
		      strc=strcmp(Name(Brothers(p)[j]),Name(pq));
		      if (strc == 0) flag[j]=1;
		      else flag[j]=0;
		      for (i=0;i<Nbrothers(pq);i++) {
			strc=strcmp(Name(Brothers(p)[j]),Name(Brothers(pq)[i]));
			if (strc == 0) flag[j]=1;
			else flag[j]=0;
		      }
		    }
		    j=0;
		    for (i=0;i<Nbrothers(p);i++) {
		      if (flag[i] == 1) {
			Brothers(p)[j]=(nodeptr) init[i];
			j++;
		      }
		    }
		    Nbrothers(p)=j;
		    twofreq(p);
		  }
		} else if (Nbrothers(p) > Nbrothers(pq)) {
		  for (i=0;i<Nbrothers(p);i++) init[i]=(bodyptr) Brothers(p)[i];
		  for (i=0;i<Nbrothers(p);i++) {
		    strc=strcmp(Name(pq),Name(Brothers(p)[i]));
		    if (strc == 0) flag[i]=0;
		    else flag[i]=1;
		  }
		  j=0;
		  for (i=0;i<Nbrothers(p);i++) {
		    if (flag[i] == 1) {
		      Brothers(p)[j]=(nodeptr) init[i];
		      j++;
		    }
		  }
		  Nbrothers(p)=j;
		  twofreq(p);

		  for (i=0;i<Nbrothers(pq);i++) init[i]=(bodyptr) Brothers(pq)[i];
		  nn1=0;
		  strc=strcmp(Name(pq),Name(p));
		  if (strc == 0) nn1=1;
		  for (i=0;i<Nbrothers(p);i++) {
		    strc=strcmp(Name(pq),Name(Brothers(p)[i]));
		    if (strc == 0) nn1=1;
		  }
		  if (nn1 == 1) {
		    for (i=0;i<Nbrothers(pq);i++) if (Frequency(pq) == Frequency(Brothers(pq)[i])) flag[i]=1;
		    else flag[i]=0;
		    j=0;
		    for (i=0;i<Nbrothers(pq);i++) {
		      if (flag[i] == 1) {
			Brothers(pq)[j]=(nodeptr) init[i];
			j++;
		      }
		    }
		    Nbrothers(pq)=j;
		    Specindx(pq)=0.0;
		    Absz(p)=0.0;
		  } else {
		    for (i=0;i<Nbrothers(pq);i++) init[i]=(bodyptr) Brothers(pq)[i];
		    for (j=0;j<Nbrothers(pq);j++) {
		      strc=strcmp(Name(Brothers(pq)[j]),Name(p));
		      if (strc == 0) flag[j]=1;
		      else flag[j]=0;
		      for (i=0;i<Nbrothers(p);i++) {
			strc=strcmp(Name(Brothers(pq)[j]),Name(Brothers(p)[i]));
			if (strc == 0) flag[j]=1;
			else flag[j]=0;
		      }
		    }
		    j=0;
		    for (i=0;i<Nbrothers(pq);i++) {
		      if (flag[i] == 1) {
			Brothers(pq)[j]=(nodeptr) init[i];
			j++;
		      }
		    }
		    Nbrothers(pq)=j;
		    twofreq(pq);
		  }
		} else if (Nbrothers(p) == Nbrothers(pq)) {
		  if (Specindx(p) > Specindx(pq)) {
		    for (i=0;i<Nbrothers(pq);i++) init[i]=(bodyptr) Brothers(pq)[i];
		    for (i=0;i<Nbrothers(pq);i++) {
		      strc=strcmp(Name(p),Name(Brothers(pq)[i]));
		      if (strc == 0) flag[i]=0;
		      else flag[i]=1;
		    }
		    j=0;
		    for (i=0;i<Nbrothers(pq);i++) {
		      if (flag[i] == 1) {
			Brothers(pq)[j]=(nodeptr) init[i];
			j++;
		      }
		    }
		    Nbrothers(pq)=j;
		    twofreq(pq);

		    for (i=0;i<Nbrothers(p);i++) init[i]=(bodyptr) Brothers(p)[i];
		    nn1=0;
		    strc=strcmp(Name(p),Name(pq));
		    if (strc == 0) nn1=1;
		    for (i=0;i<Nbrothers(pq);i++) {
		      strc=strcmp(Name(p),Name(Brothers(pq)[i]));
		      if (strc == 0) nn1=1;
		    }
		    if (nn1 == 1) {
		      for (i=0;i<Nbrothers(p);i++) if (Frequency(p) == Frequency(Brothers(p)[i])) flag[i]=1;
		      else flag[i]=0;
		      j=0;
		      for (i=0;i<Nbrothers(p);i++) {
			if (flag[i] == 1) {
			  Brothers(p)[j]=(nodeptr) init[i];
			  j++;
			}
		      }
		      Nbrothers(p)=j;
		      Specindx(p)=0.0;
		      Absz(p)=0.0;
		    } else {
		      for (i=0;i<Nbrothers(p);i++) init[i]=(bodyptr) Brothers(p)[i];
		      for (j=0;j<Nbrothers(p);j++) {
			strc=strcmp(Name(Brothers(p)[j]),Name(pq));
			if (strc == 0) flag[j]=1;
			else flag[j]=0;
			for (i=0;i<Nbrothers(pq);i++) {
			  strc=strcmp(Name(Brothers(p)[j]),Name(Brothers(pq)[i]));
			  if (strc == 0) flag[j]=1;
			  else flag[j]=0;
			}
		      }
		      j=0;
		      for (i=0;i<Nbrothers(p);i++) {
			if (flag[i] == 1) {
			  Brothers(p)[j]=(nodeptr) init[i];
			  j++;
			}
		      }
		      Nbrothers(p)=j;
		      twofreq(p);
		    }
		  } else if (Specindx(p) < Specindx(pq)) {
		    for (i=0;i<Nbrothers(p);i++) init[i]=(bodyptr) Brothers(p)[i];
		    for (i=0;i<Nbrothers(p);i++) {
		      strc=strcmp(Name(pq),Name(Brothers(p)[i]));
		      if (strc == 0) flag[i]=0;
		      else flag[i]=1;
		    }
		    j=0;
		    for (i=0;i<Nbrothers(p);i++) {
		      if (flag[i] == 1) {
			Brothers(p)[j]=(nodeptr) init[i];
			j++;
		      }
		    }
		    Nbrothers(p)=j;
		    twofreq(p);

		    for (i=0;i<Nbrothers(pq);i++) init[i]=(bodyptr) Brothers(pq)[i];
		    nn1=0;
		    strc=strcmp(Name(pq),Name(p));
		    if (strc == 0) nn1=1;
		    for (i=0;i<Nbrothers(p);i++) {
		      strc=strcmp(Name(pq),Name(Brothers(p)[i]));
		      if (strc == 0) nn1=1;
		    }
		    if (nn1 == 1) {
		      for (i=0;i<Nbrothers(pq);i++) if (Frequency(pq) == Frequency(Brothers(pq)[i])) flag[i]=1;
		      else flag[i]=0;
		      j=0;
		      for (i=0;i<Nbrothers(pq);i++) {
			if (flag[i] == 1) {
			  Brothers(pq)[j]=(nodeptr) init[i];
			  j++;
			}
		      }
		      Nbrothers(pq)=j;
		      Specindx(pq)=0.0;
		      Absz(p)=0.0;
		    } else {
		      for (i=0;i<Nbrothers(pq);i++) init[i]=(bodyptr) Brothers(pq)[i];
		      for (j=0;j<Nbrothers(pq);j++) {
			strc=strcmp(Name(Brothers(pq)[j]),Name(p));
			if (strc == 0) flag[j]=1;
			else flag[j]=0;
			for (i=0;i<Nbrothers(p);i++) {
			  strc=strcmp(Name(Brothers(pq)[j]),Name(Brothers(p)[i]));
			  if (strc == 0) flag[j]=1;
			  else flag[j]=0;
			}
		      }
		      j=0;
		      for (i=0;i<Nbrothers(pq);i++) {
			if (flag[i] == 1) {
			  Brothers(pq)[j]=(nodeptr) init[i];
			  j++;
			}
		      }
		      Nbrothers(pq)=j;
		      twofreq(pq);
		    }
		  }
		}
	      }
	    }
	  }
	}
      }
    }
  }

  for (p=bodytab;p<bodytab+nbody;p++) {
    for (ii=0;ii<Nbrothers(p);ii++) {
      pq=(bodyptr) Brothers(p)[ii];
      i=0;
	strc=strcmp(Name(pq),Name(p));
	if (strc == 0) i++;
	for (k=0;k<Nbrothers(p);k++) {
	  strc=strcmp(Name(pq),Name(Brothers(p)[k]));
	  if (strc == 0) i++;
	}
	for (j=0;j<Nbrothers(pq);j++) {
	  strc=strcmp(Name(Brothers(pq)[j]),Name(p));
	  if (strc == 0) i++;
	  for (k=0;k<Nbrothers(p);k++) {
	    strc=strcmp(Name(Brothers(pq)[j]),Name(Brothers(p)[k]));
	    if (strc == 0) i++;
	  }
	}
	if (i < Nbrothers(p)+1 || i < Nbrothers(pq)+1 || Nbrothers(p) != Nbrothers(pq)) {
	  if (i < 3 && Specindx(p)*Specindx(pq) != 0.0 && Nbrothers(p) != Nbrothers(pq)) {
	    if (Nbrothers(p) > Nbrothers(pq)) {
	      for (i=0;i<Nbrothers(pq);i++) init[i]=(bodyptr) Brothers(pq)[i];
	      for (j=0;j<Nbrothers(pq);j++) if (Frequency(Brothers(pq)[j]) == Frequency(pq)) flag[j]=1;
	      else flag[j]=0;
	      j=0;
	      for (i=0;i<Nbrothers(pq);i++) {
		if (flag[i] == 1) {
		  Brothers(pq)[j]=(nodeptr) init[i];
		  j++;
		}
	      }
	      Nbrothers(pq)=j;
	      Specindx(pq)=0.0;
	      Absz(pq)=0.0;
	      twofreq(pq);
	    } else {
	      for (i=0;i<Nbrothers(p);i++) init[i]=(bodyptr) Brothers(p)[i];
	      for (j=0;j<Nbrothers(p);j++) if (Frequency(Brothers(p)[j]) == Frequency(p)) flag[j]=1;
	      else flag[j]=0;
	      j=0;
	      for (i=0;i<Nbrothers(p);i++) {
		if (flag[i] == 1) {
		  Brothers(p)[j]=(nodeptr) init[i];
		  j++;
		}
	      }
	      Nbrothers(p)=j;
	      Specindx(p)=0.0;
	      Absz(p)=0.0;
	      twofreq(p);
	    } 
	  } else {
	    nn1=0;
	    strc=strcmp(Name(p),Name(pq));
	    if (strc == 0) nn1++;
	    for (i=0;i<Nbrothers(pq);i++) {
	      strc=strcmp(Name(p),Name(Brothers(pq)[i]));
	      if (strc == 0) nn1++;
	    }
	    if (nn1 == 0) {
	      Brothers(pq)[Nbrothers(pq)]=(nodeptr) p;
	      Nbrothers(pq)++;
	      if (Specindx(pq) == 0.0) {
		Specindx(pq)=Specindx(p);
		Absz(pq)=Absz(p);
	      }
	    }
	    nn1=0;
	    for (j=0;j<Nbrothers(p);j++) {
	      nn1=0;
	      strc=strcmp(Name(Brothers(p)[j]),Name(pq));
	      if (strc == 0) nn1++;
	      for (i=0;i<Nbrothers(pq);i++) {
		strc=strcmp(Name(Brothers(p)[j]),Name(Brothers(pq)[i]));
		if (strc == 0) nn1++;
		
	      }
	      if (nn1 == 0) {
		Brothers(pq)[Nbrothers(pq)]=Brothers(p)[j];
		Nbrothers(pq)++;
		if (Specindx(pq) == 0.0) {
		  Specindx(pq)=Specindx(p);
		  Absz(pq)=Absz(p);
		}
	      }
	    }
	    
	    nn1=0;
	    strc=strcmp(Name(pq),Name(p));
	    if (strc == 0) nn1++;
	    for (i=0;i<Nbrothers(p);i++) {
	      strc=strcmp(Name(pq),Name(Brothers(p)[i]));
	      if (strc == 0) nn1++;
	    }
	    if (nn1 == 0) {
	      Brothers(p)[Nbrothers(p)]=(nodeptr) pq;
	      Nbrothers(p)++;
	      if (Specindx(p) == 0.0) {
		Specindx(p)=Specindx(pq);
		Absz(p)=Absz(pq);
	      }
	    }
	    nn1=0;
	    for (j=0;j<Nbrothers(pq);j++) {
	      nn1=0;
	      strc=strcmp(Name(Brothers(pq)[j]),Name(p));
	      if (strc == 0) nn1++;
	      for (i=0;i<Nbrothers(p);i++) {
		strc=strcmp(Name(Brothers(pq)[j]),Name(Brothers(p)[i]));
		if (strc == 0) nn1++;
	      }
	      if (nn1 == 0) {
		Brothers(p)[Nbrothers(p)]=Brothers(pq)[j];
		Nbrothers(p)++;
		if (Specindx(p) == 0.0) {
		  Specindx(p)=Specindx(pq);
		  Absz(p)=Absz(pq);
		}
	      }
	    }
	  }
	}
     
    }
  }
  
  for (p=bodytab;p<bodytab+nbody;p++) { /* if spectrum available, but Specindx(p)=0 */
    if (Specindx(p) == 0.0 && Nbrothers(p) > 2) {
      nn1=0;
      for (ii=0;ii<Nbrothers(p);ii++) {
	pq=(bodyptr) Brothers(p)[ii];
	if (Specindx(pq) != 0.0) nn1=1; 
	if (nn1 == 1) {
	  if (Specindx(p) == 0.0) {
	    Specindx(p)=Specindx(pq);
	    Absz(p)=Absz(pq);
	  }
	  else if (Specindx(pq) < Specindx(p)) {
	    Specindx(p)=Specindx(pq);
	    Absz(p)=Absz(pq);
	  }
	}
      }
    }
  }
  
  
  for (p=bodytab;p<bodytab+nbody;p++) { /* check if Specindx != 0 and Absz = 0, then replace Absz */
    if (Specindx(p) != 0.0 && Absz(p) == 0.0) {
      for (i=0;i<Nbrothers(p);i++) {
	pq=(bodyptr) Brothers(p)[i];
	if (Specindx(pq) == Specindx(p) && Absz(pq) != 0.0) Absz(p)=Absz(pq);
      }
      if (Absz(p) == 0.0) {
	helpp=1.e20;
	for (i=0;i<Nbrothers(p);i++) {
	  pq=(bodyptr) Brothers(p)[i];
	  if (Specindx(pq) < helpp && Specindx(pq) != 0.0 && Absz(pq) != 0.0) {
	    helpp=Specindx(pq);
	    help=Absz(pq);
	  }
	}
	if (helpp != 0.0 && help != 0.0) {
	  Specindx(p)=helpp;
	  Absz(p)=help;
	}
      }
      if (Absz(p) == 0.0) {
	for (i=0;i<Nbrothers(p);i++) {
	  pq=(bodyptr) Brothers(p)[i];
	  if (Specindx(pq) != 0.0 && Absz(pq) != 0.0) {
	    Specindx(p)=Specindx(pq);
	    Absz(p)=Absz(pq);
	  }
	}
      }
    }
  }


  for (p=bodytab;p<bodytab+nbody;p++) { /* remove redundant sources */
    for (i=0;i<Nbrothers(p);i++) {
      init[i]=(bodyptr) Brothers(p)[i];
      flag[i]=1;
    }
    for (i=0;i<Nbrothers(p);i++) {
      strc=strcmp(Name(p),Name(Brothers(p)[i]));
      if (strc == 0) flag[i]=0;
    }
    for (i=0;i<Nbrothers(p)-1;i++) {
      for (j=i+1;j<Nbrothers(p);j++) {
	strc=strcmp(Name(Brothers(p)[j]),Name(Brothers(p)[i]));
	if (strc == 0) flag[i]=0;
      }
    }
    j=0;
    for (i=0;i<Nbrothers(p);i++) {
      if (flag[i] == 1) {
	Brothers(p)[j]=(nodeptr) init[i];
	j++;
      }
    }
    Nbrothers(p)=j;
  }

  for (p=bodytab;p<bodytab+nbody;p++) twofreq(p); /* security: remove source with only two different frequencies */

  
}