/*********************************************************************
 *                                                                   *
 * SPECCHECK.C  -  Checks spectra                                    *
 *                                                                   *
 * (c) 2004 by Bernd Vollmer, CDS, Observatoire de Strasbourg        *
 *                                                                   *
 *********************************************************************/
#include <stdio.h> 
#include <math.h>
#include "code.h"


void speccheck(void)
{
  int i,j, k, ik, ij, il, im, in, l, pl, ps, pb, strc, strc1, npoints, maxspec=1000, sum, specrun, pll, flag[50];
  real error, sum1, error1, factor, freq, helpp, helpq, helph, help1, help2, drab1, mor3, phii, drab,x, y, a, b, abdev, aold, bold,acalc, bcalc, freqs[50], dx, resold, dist0;
  bodyptr bbb[maxbrothers], p, pq, pp, pz, p1, px, pq1, pq2;
  bodyptr spectrum[maxspec], init[51];
  real nu, flux, fluxerror, logflux, lognu, estflux, spec[maxspec];
  vector pos0;
 
 
  factor=1.5;                  /* 1.5 fudge factor for flux error */

  for (p=bodytab;p<bodytab+nbody;p++) {

    factor=1.5;
    specrun=1;
    if (Specindx(p) != 0.0) {

    again:

      pl=Nbrothers(p);
      k=0;                 /* keep brothers at same frequency */
      freq=Frequency(p);
      if (Nbrothers(p) > 0) {
	for (i=0;i<Nbrothers(p);i++) bbb[i]=(bodyptr) Brothers(p)[i];
	for (i=0;i<Nbrothers(p);i++) {
	  if (Frequency(bbb[i]) == freq) {
	    Brothers(p)[k]=(nodeptr) bbb[i];
	    k++;
	  }
	}
	Nbrothers(p)=k;
		
	
	j=0;
	for (i=0;i<pl;i++) {
	  nu=Frequency(Brothers(p)[i]);
	  if (nu != freq) {                            /* at different frequency */
	    helpp=Absz(p)+Specindx(p)*log10(nu); 
	    estflux=pow(10.,helpp);
	    flux=Flux(Brothers(p)[i]);
	    fluxerror=Fluxerror(Brothers(p)[i]);
	    if (estflux <= flux+factor*fluxerror &&  estflux >= flux-factor*fluxerror) {
	      spectrum[j]=(bodyptr) Brothers(p)[i];
	      j++;
	    }
	  }
	}
	
	npoints=j;
	
	pll=0;                        /* apply proximity criterion */
	pp=p;
	for (i=0;i<npoints;i++) { 
	  pq=spectrum[i];
	  if (Majaxis(pp) == 0) helpp=Resolution(pp);
	  else helpp=Majaxis(pp);
	  if (Majaxis(pq) == 0) helpq=Resolution(pq);
	  else helpq=Majaxis(pq);
	  helpp+=helpq;
	  helpp/=2.;
	  DISTV1(helph,Pos(pp),Pos(pq));
	  if (helph > helpp) pll=1;
	}
	for (i=0;i<npoints-1;i++) { 
	  for (j=i+1;j<npoints;j++) {
	    pq=spectrum[i];
	    pp=spectrum[j];
	    if (Majaxis(pp) == 0) helpp=Resolution(pp);
	    else helpp=Majaxis(pp);
	    if (Majaxis(pq) == 0) helpq=Resolution(pq);
	    else helpq=Majaxis(pq);
	    helpp+=helpq;
	    helpp/=2.;
	    DISTV1(helph,Pos(pp),Pos(pq));
	    if (helph > helpp) pll=1;
	  }
	}
	if (pll == 1) {
	  Specindx(p)=0.0;
	  Absz(p)=0.0;
	} else {
	  
	  pb=0;                       /* initial source must fit into the spectrum */
	  nu=Frequency(p);
	  helpp=Absz(p)+Specindx(p)*log10(nu); 
	  estflux=pow(10.,helpp);
	  flux=Flux(p);
	  fluxerror=Fluxerror(p);
	  if (estflux <= flux+factor*fluxerror &&  estflux >= flux-factor*fluxerror) pb=1;
	  
	  if (npoints > 1 && pb == 1) {
	    for (i=0;i<npoints;i++) {
	      j=Nbrothers(p)+i;
	      Brothers(p)[j]=(nodeptr) spectrum[i];
	    }
	    Nbrothers(p)+=npoints;
	  } 
	  
	  j=0;                  /* count number of different frequencies */
	  for (i=0;i<Nbrothers(p);i++) {
	    nu=Frequency(Brothers(p)[i]);
	    if (nu != freq) {
	      spec[j]=Frequency(Brothers(p)[i]);
	      j++;
	    }
	  }
	  for (i=0;i<j;i++) {
	    for (k=0;k<j;k++) {
	      if (i != k && spec[i] == spec[k]) spec[i]=0.0;
	    }
	  }
	  sum=0;
	  for (i=0;i<j;i++) if (spec[i] != 0.0) sum++;   /* sum > 1 for retaining of spectrum */
	  if (sum < 2) {
	    Specindx(p)=0.0;
	    Absz(p)=0.0;
	  } else {
	    
	    spec[j]=Frequency(p);
	    j++;
	    for (i=0;i<j;i++) freqs[i]=1.e21;
	    for (k=0;k<j;k++) {
	      l=-1;
	      for (i=0;i<j;i++) {
		if (spec[i] != 0.0 && spec[i] < freqs[k]) {
		  freqs[k]=spec[i];
		  l=i;
		}
	      }
	      if (l >= 0) spec[l]=0.0;
	    }
	    l=0;
	    for (i=0;i<j;i++) if (freqs[i] < 1.e20) l++;
	    for (i=1;i<l;i++) spec[i]=freqs[i]-freqs[i-1];
	    sum1=0.;
	    k=-1;
	    for (i=1;i<l;i++) {
	      if (spec[i] > sum1) {
		sum1=spec[i];
		k=i;
	      }
	    }
	    if (k >= 0) spec[k]=0.0;
	    else {
	      printf("problem in speccheck \n");
	      getchar();
	    }
	    error1=0.0;
	    for (i=1;i<l;i++) error1=MAX(error1,spec[i]);
	    if (error1 != 0.0) helpp=error1/sum1;
	    if (sum1 <= 0.0 || error1 <= 0.0 || helpp > 1.) {
	      printf("ratio=%5.3e  error in speccheck\n",helpp);
	      getchar();
	    }
	    if (specrun == 1) {
	      if (helpp < 0.1) {
		factor=1.0;
		specrun++;
		goto again;
	      }
	    }
	    if (helpp < 0.02) {
	      Specindx(p)=0.0;
	      Absz(p)=0.0;
	    }
	  }
	}
      }
    }
  }

  printf("beginning long check...\n");
  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");
  
  pl=ps=0;
  for (p=bodytab;p<bodytab+nbody;p++) {
    resold=Resolution(p);
    Resolution(p)=1200.;
    hackgrav1(p, Mass(p) > 0.0);    /* proximity search */	  
    Resolution(p)=resold;
    /*printf("ccount=%i \n",ccount);*/
    ps++;
    if (Nbrothers(p) > 0 && ccount > 0) {
      for (i=0;i<Nbrothers(p);i++) {
	pq=(bodyptr) Brothers(p)[i];
	for (in=0;in<ccount;in++) {
	  p1=(bodyptr) store[in];
	  for (ij=0;ij<Nbrothers(p1);ij++) {
	    pq1=(bodyptr) Brothers(p1)[ij];
	    strc=strcmp(Name(pq),Name(pq1));
	    if (strc == 0) {
	      strc1=1;
	      for (il=0;il<Nbrothers(p);il++) if (il != i) strc1*=strncmp(Name(Brothers(p)[il]),Name(p1),3);
	      for (im=0;im<Nbrothers(p1);im++) if (im != ij) strc1*=strncmp(Name(p),Name(Brothers(p1)[im]),3);
	      for (il=0;il<Nbrothers(p);il++) {
		for (im=0;im<Nbrothers(p1);im++) {
		  if (il != i && im != ij) strc1*=strncmp(Name(Brothers(p)[il]),Name(Brothers(p1)[im]),3);
		}
	      }
	      if (strc1 == 0) {
		strc1=1;
		for (il=0;il<Nbrothers(p);il++) if (il != i) strc1*=strcmp(Name(Brothers(p)[il]),Name(p1));
		for (im=0;im<Nbrothers(p1);im++) if (im != ij) strc1*=strcmp(Name(p),Name(Brothers(p1)[im]));
		for (il=0;il<Nbrothers(p);il++) {
		  for (im=0;im<Nbrothers(p1);im++) {
		    if (il != i && im != ij) strc1*=strcmp(Name(Brothers(p)[il]),Name(Brothers(p1)[im]));
		  }
		}
		if (strc1 != 0) {
		  help1=1./Resolution(p)*Pos(p)[0];
		  help2=1./Resolution(p)*Pos(p)[1];
		  helph=1./Resolution(p);
		  for (ik=0;ik<Nbrothers(p);ik++) {
		    pq2=(bodyptr) Brothers(p)[ik];
		    help1+=1./Resolution(pq2)*Pos(pq2)[0];
		    help2+=1./Resolution(pq2)*Pos(pq2)[1];
		    helph+=1./Resolution(pq2);
		  }
		  help1/=helph;
		  help2/=helph;
		  pos0[0]=help1;
		  pos0[1]=help2;
		  DISTV1(helpp,pos0,Pos(pq));
		  help1=1./Resolution(p1)*Pos(p1)[0];
		  help2=1./Resolution(p1)*Pos(p1)[1];
		  helph=1./Resolution(p1);
		  for (ik=0;ik<Nbrothers(p1);ik++) {
		    pq2=(bodyptr) Brothers(p1)[ik];
		    help1+=1./Resolution(pq2)*Pos(pq2)[0];
		    help2+=1./Resolution(pq2)*Pos(pq2)[1];
		    helph+=1./Resolution(pq2);
		  }
		  help1/=helph;
		  help2/=helph;
		  pos0[0]=help1;
		  pos0[1]=help2;
		  DISTV1(helpq,pos0,Pos(pq1));
		  if (helpp < helpq) {
		    for (k=0;k<Nbrothers(p1);k++) bbb[k]=(bodyptr) Brothers(p1)[k];
		    l=0;
		    for (k=0;k<Nbrothers(p1);k++) {
		      if (k != ij) {
			Brothers(p1)[l]=(nodeptr) bbb[k];
			l++;
		      }
		    }
		    pl++;
		    if (l != Nbrothers(p1)-1) {
		      printf("problem in speccheck %s %s case1 \n",Name(pq),Name(pq1));
		    }
		    Nbrothers(p1)=l;
		  }
		  if (helpp > helpq) {
		    for (k=0;k<Nbrothers(p);k++) bbb[k]=(bodyptr) Brothers(p)[k];
		    l=0;
		    for (k=0;k<Nbrothers(p);k++) {
		      if (k != i) {
			Brothers(p)[l]=(nodeptr) bbb[k];
			l++;
		      }
		    }
		    pl++;
		    Nbrothers(p)=l;
		  }
		}
	      }
	    }
	  }
	}
      }
    }
  }

  if (doedge) {
    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;
      }
    }
    for (p=bodytab;p<bodytab+nbody;p++) {
      if (Mass(p) == 2.) Pos(p)[0]+=dist0;
      if (Mass(p) == 3.) Pos(p)[0]-=dist0;
      resold=Resolution(p);
      Resolution(p)=1200.;
      hackgrav1(p, Mass(p) > 0.0);    /* proximity search */	  
      Resolution(p)=resold;
      if (Mass(p) == 2.) Pos(p)[0]-=dist0;
      if (Mass(p) == 3.) Pos(p)[0]+=dist0;
      ps++;
      if (Nbrothers(p) > 0 && ccount > 0) {
	for (i=0;i<Nbrothers(p);i++) {
	  pq=(bodyptr) Brothers(p)[i];
	  for (in=0;in<ccount;in++) {
	    p1=(bodyptr) store[in];
	    for (ij=0;ij<Nbrothers(p1);ij++) {
	      pq1=(bodyptr) Brothers(p1)[ij];
	      strc=strcmp(Name(pq),Name(pq1));
	      if (strc == 0) {
		strc1=1;
		for (il=0;il<Nbrothers(p);il++) if (il != i) strc1*=strncmp(Name(Brothers(p)[il]),Name(p1),3);
		for (im=0;im<Nbrothers(p1);im++) if (im != ij) strc1*=strncmp(Name(p),Name(Brothers(p1)[im]),3);
		for (il=0;il<Nbrothers(p);il++) {
		  for (im=0;im<Nbrothers(p1);im++) {
		    if (il != i && im != ij) strc1*=strncmp(Name(Brothers(p)[il]),Name(Brothers(p1)[im]),3);
		  }
		}
		if (strc1 == 0) {
		  strc1=1;
		  for (il=0;il<Nbrothers(p);il++) if (il != i) strc1*=strcmp(Name(Brothers(p)[il]),Name(p1));
		  for (im=0;im<Nbrothers(p1);im++) if (im != ij) strc1*=strcmp(Name(p),Name(Brothers(p1)[im]));
		  for (il=0;il<Nbrothers(p);il++) {
		    for (im=0;im<Nbrothers(p1);im++) {
		      if (il != i && im != ij) strc1*=strcmp(Name(Brothers(p)[il]),Name(Brothers(p1)[im]));
		    }
		  }
		  if (strc1 != 0) {
		    help1=1./Resolution(p)*Pos(p)[0];
		    help2=1./Resolution(p)*Pos(p)[1];
		    helph=1./Resolution(p);
		    for (ik=0;ik<Nbrothers(p);ik++) {
		      pq2=(bodyptr) Brothers(p)[ik];
		      help1+=1./Resolution(pq2)*Pos(pq2)[0];
		      help2+=1./Resolution(pq2)*Pos(pq2)[1];
		      helph+=1./Resolution(pq2);
		    }
		    help1/=helph;
		    help2/=helph;
		    pos0[0]=help1;
		    pos0[1]=help2;
		    DISTV1(helpp,pos0,Pos(pq));
		    help1=1./Resolution(p1)*Pos(p1)[0];
		    help2=1./Resolution(p1)*Pos(p1)[1];
		    helph=1./Resolution(p1);
		    for (ik=0;ik<Nbrothers(p1);ik++) {
		      pq2=(bodyptr) Brothers(p1)[ik];
		      help1+=1./Resolution(pq2)*Pos(pq2)[0];
		      help2+=1./Resolution(pq2)*Pos(pq2)[1];
		      helph+=1./Resolution(pq2);
		    }
		    help1/=helph;
		    help2/=helph;
		    pos0[0]=help1;
		    pos0[1]=help2;
		    DISTV1(helpq,pos0,Pos(pq1));
		    if (helpp < helpq) {
		      for (k=0;k<Nbrothers(p1);k++) bbb[k]=(bodyptr) Brothers(p1)[k];
		      l=0;
		      for (k=0;k<Nbrothers(p1);k++) {
			if (k != ij) {
			  Brothers(p1)[l]=(nodeptr) bbb[k];
			  l++;
			}
		      }
		      pl++;
		      if (l != Nbrothers(p1)-1) {
			printf("problem in speccheck %s %s case1 \n",Name(pq),Name(pq1));
		      }
		      Nbrothers(p1)=l;
		    }
		    if (helpp > helpq) {
		      for (k=0;k<Nbrothers(p);k++) bbb[k]=(bodyptr) Brothers(p)[k];
		      l=0;
		      for (k=0;k<Nbrothers(p);k++) {
			if (k != i) {
			  Brothers(p)[l]=(nodeptr) bbb[k];
			  l++;
			}
		      }
		      pl++;
		      Nbrothers(p)=l;
		    }
		  }
		}
	      }
	    }
	  }
	}
      }
    }
    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;
    }
  }
  printf("%i ambiguous sources removed \n",pl);


  for (p=bodytab;p<bodytab+nbody;p++) {
    if (Specindx(p) == 0.0 && Absz(p) == 0.0) {
      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)) 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;
    }
  }

}