/*********************************************************************
 *                                                                   *
 * SPECTRA.C  -  Extracts spectra                                    *
 *                                                                   *
 * (c) 2004 by Bernd Vollmer, CDS, Observatoire de Strasbourg        *
 *                                                                   *
 *********************************************************************/
#include <stdio.h>
#include "nrutil.h" 
#include "nrutil.c" 
#include <math.h>
#include "code.h"
#include "medfit.c"
#include "fit.c"

local int ffcompare(bodyptr,bodyptr);
extern real em1, em2;

void spectra(bodyptr p)
{
  FILE *datei_ptr;
  int i,j, k, pl, ps, pb, strc, npoints, npoints1, npoints2, npoints3, *spflag, *wat, maxspec=1000, qwe, kk, khelp, pll;
  short check, check1, check2, check3, check4, what[maxspec], what1[maxspec], zaehler, zaehler1;
  real sum, error, sum1, error1, factor, freq, helpp, helpq, helph, drab1, mor3, phii, drab,x, y, a, b, abdev, aold, bold,acalc, bcalc, chi2, siga, sigb, freqhelp[51], dx;
  bodyptr ppp[maxparents],sss[maxsons],bbb[maxbrothers], pq;
  short pflag[maxparents],bflag[maxbrothers],sflag[maxsons];
  bodyptr spectrum[maxspec], spectrum1[maxspec], spectrum2[maxspec], pskip, q, pp;
  real *nu, *flux, *fluxerror, *logflux, *lognu, *estflux;
  char Datei[11], ccomp[20], asd[20];
  comppareptr qq;
  nu=vector0(1,maxspec);
  flux=vector0(1,maxspec);
  fluxerror=vector0(1,maxspec);
  logflux=vector0(1,maxspec);
  lognu=vector0(1,maxspec);
  estflux=vector0(1,maxspec);
  spflag=ivector0(1,maxspec);
  wat=ivector0(1,maxspec);

  zaehler=zaehler1=0;
  factor=1.5;                  /* 1.5 fudge factor for flux error */

  /* begin of dependences restoring */

  j=0;
  what[j]=SOURCE;              /* source description */
  spectrum[j]=p;               /* point in spectrum */
  j++;

  freq=Frequency(p);
  for (i=0;i<Nbrothers(p);i++) {      /* restore brother dependences */
    if (Frequency(Brothers(p)[i]) != freq) { /* at different frequency */
      what[j]=BROTHER;
      pq=(bodyptr) Brothers(p)[i];    /* check brothers */
      spectrum[j]=pq;
      j++;
      pskip=p;
      for (k=0;k<Nparents(pq);k++) {
	q=(bodyptr) Parents(pq)[k];   /* parents of brother */
	if (Frequency(pq) == Frequency(q)) {  /* at the same frequency */
	  check=check1=0;
	  helpp=MAX(Resolution(pskip),Majaxis(pskip));
	  helpq=MAX(Resolution(q),Majaxis(q));
	  
	  if (Majaxis(pskip) == 0) helpp=Resolution(pskip);
	  if (Majaxis(q) == 0) helpq=Resolution(q);
	  
	jump:
	  
	  if (helpq > helpp*1.25) {
	    if (Resolution(q) >= Majaxis(q)) {
	      DISTV(helph,Pos(pskip),Pos(q));
	      drab1=Resolution(q);
	      if (helph < drab1/2.) {
		spectrum[j]=q;
		what[j]=POB;          /* parent of brother */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    } else {
	      x=Pos(pskip)[0]-Pos(q)[0];
	      y=Pos(pskip)[1]-Pos(q)[1];
	      phii=(PA(q))/180.*PI;
	      mor3=(x*cos(phii)-y*sin(phii))*(x*cos(phii)-y*sin(phii))/(Majaxis(q)/2.*Majaxis(q)/2.)+(x*sin(phii)+y*cos(phii))*(x*sin(phii)+y*cos(phii))/(Minaxis(q)/2.*Minaxis(q)/2.);
	      if (mor3 < 1.) {
		spectrum[j]=q;
		what[j]=POB;          /* parent of brother */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    }	
	  }
	  else if (helpq < 0.75*helpp) { 
	    if (Resolution(pskip) >= Majaxis(pskip)) {
	      DISTV(helph,Pos(pskip),Pos(q));
	      drab1=Resolution(pskip);
	      if (helph < drab1/2.) {
		spectrum[j]=q;
		what[j]=SOB;           /* son of brother */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    } else {
	      x=Pos(q)[0]-Pos(pskip)[0];
	      y=Pos(q)[1]-Pos(pskip)[1];
	      phii=(PA(pskip))/180.*PI;
	      mor3=(x*cos(phii)-y*sin(phii))*(x*cos(phii)-y*sin(phii))/(Majaxis(pskip)/2.*Majaxis(pskip)/2.)+(x*sin(phii)+y*cos(phii))*(x*sin(phii)+y*cos(phii))/(Minaxis(pskip)/2.*Minaxis(pskip)/2.);
	      if (mor3 < 1.) {
		spectrum[j]=q;
		what[j]=SOB;            /* son of brother */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    }	
	  }
	  drab=MAX(helpp,helpq);

	  /* begin check around RA = 00 00 00 */

	  if (Pos(pskip)[0]+2.0*drab > 1.296e6 && check == 0 && check1 == 0) { 
	    Pos(q)[0]-=1.296e6;
	    check=1;
	    goto jump;
	  }
	  if (Pos(pskip)[0]-2.0*drab < 0. && check == 0 && check1 == 0) {
	    Pos(q)[0]+=1.296e6;
	    check=2;
	    goto jump;
	  }
	  if (check == 1) Pos(q)[0]+=1.296e6;
	  if (check == 2) Pos(q)[0]-=1.296e6;  
	  
	  /* end check around RA = 00 00 00 */

	}
      }

      for (k=0;k<Nsons(pq);k++) {
	q=(bodyptr) Sons(pq)[k];   /* sons of brother */
	if (Frequency(pq) == Frequency(q)) {  /* at the same frequency */
	  check=check1=0;
	  helpp=MAX(Resolution(pskip),Majaxis(pskip));
	  helpq=MAX(Resolution(q),Majaxis(q));
	  
	  if (Majaxis(pskip) == 0) helpp=Resolution(pskip);
	  if (Majaxis(q) == 0) helpq=Resolution(q);
	  
	jumpp:
	  
	  if (helpq > helpp*1.25) {
	    if (Resolution(q) >= Majaxis(q)) {
	      DISTV(helph,Pos(pskip),Pos(q));
	      drab1=Resolution(q);
	      if (helph < drab1/2.) {
		spectrum[j]=q;
		what[j]=POB;          /* parent of brother */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    } else {
	      x=Pos(pskip)[0]-Pos(q)[0];
	      y=Pos(pskip)[1]-Pos(q)[1];
	      phii=(PA(q))/180.*PI;
	      mor3=(x*cos(phii)-y*sin(phii))*(x*cos(phii)-y*sin(phii))/(Majaxis(q)/2.*Majaxis(q)/2.)+(x*sin(phii)+y*cos(phii))*(x*sin(phii)+y*cos(phii))/(Minaxis(q)/2.*Minaxis(q)/2.);
	      if (mor3 < 1.) {
		spectrum[j]=q;
		what[j]=POB;          /* parent of brother */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    }	
	  }
	  else if (helpq < 0.75*helpp) { 
	    if (Resolution(pskip) >= Majaxis(pskip)) {
	      DISTV(helph,Pos(pskip),Pos(q));
	      drab1=Resolution(pskip);
	      if (helph < drab1/2.) {
		spectrum[j]=q;
		what[j]=SOB;           /* son of brother */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    } else {
	      x=Pos(q)[0]-Pos(pskip)[0];
	      y=Pos(q)[1]-Pos(pskip)[1];
	      phii=(PA(pskip))/180.*PI;
	      mor3=(x*cos(phii)-y*sin(phii))*(x*cos(phii)-y*sin(phii))/(Majaxis(pskip)/2.*Majaxis(pskip)/2.)+(x*sin(phii)+y*cos(phii))*(x*sin(phii)+y*cos(phii))/(Minaxis(pskip)/2.*Minaxis(pskip)/2.);
	      if (mor3 < 1.) {
		spectrum[j]=q;
		what[j]=SOB;            /* son of brother */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    }	
	  }
	  drab=MAX(helpp,helpq);

	  /* begin check around RA = 00 00 00 */

	  if (Pos(pskip)[0]+2.0*drab > 1.296e6 && check == 0 && check1 == 0) { 
	    Pos(q)[0]-=1.296e6;
	    check=1;
	    goto jumpp;
	  }
	  if (Pos(pskip)[0]-2.0*drab < 0. && check == 0 && check1 == 0) {
	    Pos(q)[0]+=1.296e6;
	    check=2;
	    goto jumpp;
	  }
	  if (check == 1) Pos(q)[0]+=1.296e6;
	  if (check == 2) Pos(q)[0]-=1.296e6;  
	  
	  /* end check around RA = 00 00 00 */

	}
      }
    }
  }
  
  for (i=0;i<Nparents(p);i++) {           /* restore parent dependences */
    if (Frequency(Parents(p)[i]) != freq) { /* at different frequency */
      what[j]=PARENT;
      pq=(bodyptr) Parents(p)[i];         /* check parents */
      spectrum[j]=pq;
      j++;
      pskip=p;
      for (k=0;k<Nparents(pq);k++) {
	q=(bodyptr) Parents(pq)[k];       /* parents of parents */
	if (Frequency(pq) == Frequency(q)) { /* at the same frequency */
	  check=check1=0;
	  helpp=MAX(Resolution(pskip),Majaxis(pskip));
	  helpq=MAX(Resolution(q),Majaxis(q));
	  
	  if (Majaxis(pskip) == 0) helpp=Resolution(pskip);
	  if (Majaxis(q) == 0) helpq=Resolution(q);
	  
	jump1:
	  
	  if (helpq > helpp*1.25) {
	    if (Resolution(q) >= Majaxis(q)) {
	      DISTV(helph,Pos(pskip),Pos(q));
	      drab1=Resolution(q);
	      if (helph < drab1/2.) {
		spectrum[j]=q;
		what[j]=POP;                /* parent of parent */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    } else {
	      x=Pos(pskip)[0]-Pos(q)[0];
	      y=Pos(pskip)[1]-Pos(q)[1];
	      phii=(PA(q))/180.*PI;
	      mor3=(x*cos(phii)-y*sin(phii))*(x*cos(phii)-y*sin(phii))/(Majaxis(q)/2.*Majaxis(q)/2.)+(x*sin(phii)+y*cos(phii))*(x*sin(phii)+y*cos(phii))/(Minaxis(q)/2.*Minaxis(q)/2.);
	      if (mor3 < 1.) {
		spectrum[j]=q;
		what[j]=POP;                 /* parent of parent */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    }	
	  }
	  else if (helpq < 0.75*helpp) { 
	    if (Resolution(pskip) >= Majaxis(pskip)) {
	      DISTV(helph,Pos(pskip),Pos(q));
	      drab1=Resolution(pskip);
	      if (helph < drab1/2.) {
		spectrum[j]=q;
		what[j]=SOP;                  /* son of parent */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    } else {
	      x=Pos(q)[0]-Pos(pskip)[0];
	      y=Pos(q)[1]-Pos(pskip)[1];
	      phii=(PA(pskip))/180.*PI;
	      mor3=(x*cos(phii)-y*sin(phii))*(x*cos(phii)-y*sin(phii))/(Majaxis(pskip)/2.*Majaxis(pskip)/2.)+(x*sin(phii)+y*cos(phii))*(x*sin(phii)+y*cos(phii))/(Minaxis(pskip)/2.*Minaxis(pskip)/2.);
	      if (mor3 < 1.) {
		spectrum[j]=q;
		what[j]=SOP;                  /* son of parent */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    }	
	  }

	  /* begin check around RA = 00 00 00 */

	  drab=MAX(helpp,helpq);
	  if (Pos(pskip)[0]+2.0*drab > 1.296e6 && check == 0 && check1 == 0) {
	    Pos(q)[0]-=1.296e6;
	    check=1;
	    goto jump1;
	  }
	  if (Pos(pskip)[0]-2.0*drab < 0. && check == 0 && check1 == 0) {
	    Pos(q)[0]+=1.296e6;
	    check=2;
	    goto jump1;
	  }
	  if (check == 1) Pos(q)[0]+=1.296e6;
	  if (check == 2) Pos(q)[0]-=1.296e6;   

	  /* end check around RA = 00 00 00 */

	}
      }

      for (k=0;k<Nsons(pq);k++) {
	q=(bodyptr) Sons(pq)[k];       /* sons of parents */
	if (Frequency(pq) == Frequency(q)) { /* at the same frequency */
	  check=check1=0;
	  helpp=MAX(Resolution(pskip),Majaxis(pskip));
	  helpq=MAX(Resolution(q),Majaxis(q));
	  
	  if (Majaxis(pskip) == 0) helpp=Resolution(pskip);
	  if (Majaxis(q) == 0) helpq=Resolution(q);
	  
	jumpp1:
	  
	  if (helpq > helpp*1.25) {
	    if (Resolution(q) >= Majaxis(q)) {
	      DISTV(helph,Pos(pskip),Pos(q));
	      drab1=Resolution(q);
	      if (helph < drab1/2.) {
		spectrum[j]=q;
		what[j]=POP;                /* parent of parent */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    } else {
	      x=Pos(pskip)[0]-Pos(q)[0];
	      y=Pos(pskip)[1]-Pos(q)[1];
	      phii=(PA(q))/180.*PI;
	      mor3=(x*cos(phii)-y*sin(phii))*(x*cos(phii)-y*sin(phii))/(Majaxis(q)/2.*Majaxis(q)/2.)+(x*sin(phii)+y*cos(phii))*(x*sin(phii)+y*cos(phii))/(Minaxis(q)/2.*Minaxis(q)/2.);
	      if (mor3 < 1.) {
		spectrum[j]=q;
		what[j]=POP;                 /* parent of parent */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    }	
	  }
	  else if (helpq < 0.75*helpp) { 
	    if (Resolution(pskip) >= Majaxis(pskip)) {
	      DISTV(helph,Pos(pskip),Pos(q));
	      drab1=Resolution(pskip);
	      if (helph < drab1/2.) {
		spectrum[j]=q;
		what[j]=SOP;                  /* son of parent */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    } else {
	      x=Pos(q)[0]-Pos(pskip)[0];
	      y=Pos(q)[1]-Pos(pskip)[1];
	      phii=(PA(pskip))/180.*PI;
	      mor3=(x*cos(phii)-y*sin(phii))*(x*cos(phii)-y*sin(phii))/(Majaxis(pskip)/2.*Majaxis(pskip)/2.)+(x*sin(phii)+y*cos(phii))*(x*sin(phii)+y*cos(phii))/(Minaxis(pskip)/2.*Minaxis(pskip)/2.);
	      if (mor3 < 1.) {
		spectrum[j]=q;
		what[j]=SOP;                  /* son of parent */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    }	
	  }

	  /* begin check around RA = 00 00 00 */

	  drab=MAX(helpp,helpq);
	  if (Pos(pskip)[0]+2.0*drab > 1.296e6 && check == 0 && check1 == 0) {
	    Pos(q)[0]-=1.296e6;
	    check=1;
	    goto jumpp1;
	  }
	  if (Pos(pskip)[0]-2.0*drab < 0. && check == 0 && check1 == 0) {
	    Pos(q)[0]+=1.296e6;
	    check=2;
	    goto jumpp1;
	  }
	  if (check == 1) Pos(q)[0]+=1.296e6;
	  if (check == 2) Pos(q)[0]-=1.296e6;   

	  /* end check around RA = 00 00 00 */

	}
      }
    }
  }
  
  for (i=0;i<Nsons(p);i++) {              /* restore son dependeces */
    if (Frequency(Sons(p)[i]) != freq) {  /* at different frequency */
      what[j]=SON;
      pq=(bodyptr) Sons(p)[i];            /* check sons */
      spectrum[j]=pq;
      j++;
      pskip=p;
      for (k=0;k<Nparents(pq);k++) {
	q=(bodyptr) Parents(pq)[k];       /* parents of sons */
	if (Frequency(pq) == Frequency(q)) { /* at the same frequency */
	  check=check1=0;
	  helpp=MAX(Resolution(pskip),Majaxis(pskip));
	  helpq=MAX(Resolution(q),Majaxis(q));
	  
	  if (Majaxis(pskip) == 0) helpp=Resolution(pskip);
	  if (Majaxis(q) == 0) helpq=Resolution(q);
	  
	jump2:
	  
	  if (helpq > helpp*1.25) {
	    if (Resolution(q) >= Majaxis(q)) {
	      DISTV(helph,Pos(pskip),Pos(q));
	      drab1=Resolution(q);
	      if (helph < drab1/2.) {
		spectrum[j]=q;
		what[j]=POS;              /* parent of son */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    } else {
	      x=Pos(pskip)[0]-Pos(q)[0];
	      y=Pos(pskip)[1]-Pos(q)[1];
	      phii=(PA(q))/180.*PI;
	      mor3=(x*cos(phii)-y*sin(phii))*(x*cos(phii)-y*sin(phii))/(Majaxis(q)/2.*Majaxis(q)/2.)+(x*sin(phii)+y*cos(phii))*(x*sin(phii)+y*cos(phii))/(Minaxis(q)/2.*Minaxis(q)/2.);
	      if (mor3 < 1.) {
		spectrum[j]=q;
		what[j]=POS;              /* parent of son */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    }	
	  }
	  else if (helpq < 0.75*helpp) { 
	    if (Resolution(pskip) >= Majaxis(pskip)) {
	      DISTV(helph,Pos(pskip),Pos(q));
	      drab1=Resolution(pskip);
	      if (helph < drab1/2.) {
		spectrum[j]=q;
		what[j]=SOS;               /* son of son */  
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    } else {
	      x=Pos(q)[0]-Pos(pskip)[0];
	      y=Pos(q)[1]-Pos(pskip)[1];
	      phii=(PA(pskip))/180.*PI;
	      mor3=(x*cos(phii)-y*sin(phii))*(x*cos(phii)-y*sin(phii))/(Majaxis(pskip)/2.*Majaxis(pskip)/2.)+(x*sin(phii)+y*cos(phii))*(x*sin(phii)+y*cos(phii))/(Minaxis(pskip)/2.*Minaxis(pskip)/2.);
	      if (mor3 < 1.) {
		spectrum[j]=q;
		what[j]=SOS;               /* son of son */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    }	
	  }

	  /* begin check around RA = 00 00 00 */

	  drab=MAX(helpp,helpq);
	  if (Pos(pskip)[0]+2.0*drab > 1.296e6 && check == 0 && check1 == 0) {
	    Pos(q)[0]-=1.296e6;
	    check=1;
	    goto jump2;
	  }
	  if (Pos(pskip)[0]-2.0*drab < 0. && check == 0 && check1 == 0) {
	    Pos(q)[0]+=1.296e6;
	    check=2;
	    goto jump2;
	  }
	  if (check == 1) Pos(q)[0]+=1.296e6;
	  if (check == 2) Pos(q)[0]-=1.296e6;

	  /* end check around RA = 00 00 00 */
   
	}
      }

      for (k=0;k<Nsons(pq);k++) {
	q=(bodyptr) Sons(pq)[k];       /* sons of sons */
	if (Frequency(pq) == Frequency(q)) { /* at the same frequency */
	  check=check1=0;
	  helpp=MAX(Resolution(pskip),Majaxis(pskip));
	  helpq=MAX(Resolution(q),Majaxis(q));
	  
	  if (Majaxis(pskip) == 0) helpp=Resolution(pskip);
	  if (Majaxis(q) == 0) helpq=Resolution(q);
	  
	jumpp2:
	  
	  if (helpq > helpp*1.25) {
	    if (Resolution(q) >= Majaxis(q)) {
	      DISTV(helph,Pos(pskip),Pos(q));
	      drab1=Resolution(q);
	      if (helph < drab1/2.) {
		spectrum[j]=q;
		what[j]=POS;              /* parent of son */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    } else {
	      x=Pos(pskip)[0]-Pos(q)[0];
	      y=Pos(pskip)[1]-Pos(q)[1];
	      phii=(PA(q))/180.*PI;
	      mor3=(x*cos(phii)-y*sin(phii))*(x*cos(phii)-y*sin(phii))/(Majaxis(q)/2.*Majaxis(q)/2.)+(x*sin(phii)+y*cos(phii))*(x*sin(phii)+y*cos(phii))/(Minaxis(q)/2.*Minaxis(q)/2.);
	      if (mor3 < 1.) {
		spectrum[j]=q;
		what[j]=POS;              /* parent of son */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    }	
	  }
	  else if (helpq < 0.75*helpp) { 
	    if (Resolution(pskip) >= Majaxis(pskip)) {
	      DISTV(helph,Pos(pskip),Pos(q));
	      drab1=Resolution(pskip);
	      if (helph < drab1/2.) {
		spectrum[j]=q;
		what[j]=SOS;               /* son of son */  
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    } else {
	      x=Pos(q)[0]-Pos(pskip)[0];
	      y=Pos(q)[1]-Pos(pskip)[1];
	      phii=(PA(pskip))/180.*PI;
	      mor3=(x*cos(phii)-y*sin(phii))*(x*cos(phii)-y*sin(phii))/(Majaxis(pskip)/2.*Majaxis(pskip)/2.)+(x*sin(phii)+y*cos(phii))*(x*sin(phii)+y*cos(phii))/(Minaxis(pskip)/2.*Minaxis(pskip)/2.);
	      if (mor3 < 1.) {
		spectrum[j]=q;
		what[j]=SOS;               /* son of son */
		j++;
		if (j >= maxspec) {
		  printf("too many spectral points \n");
		  exit(0);
		}
		check1=1;
	      }   
	    }	
	  }

	  /* begin check around RA = 00 00 00 */

	  drab=MAX(helpp,helpq);
	  if (Pos(pskip)[0]+2.0*drab > 1.296e6 && check == 0 && check1 == 0) {
	    Pos(q)[0]-=1.296e6;
	    check=1;
	    goto jumpp2;
	  }
	  if (Pos(pskip)[0]-2.0*drab < 0. && check == 0 && check1 == 0) {
	    Pos(q)[0]+=1.296e6;
	    check=2;
	    goto jumpp2;
	  }
	  if (check == 1) Pos(q)[0]+=1.296e6;
	  if (check == 2) Pos(q)[0]-=1.296e6;

	  /* end check around RA = 00 00 00 */
   
	}
      }
    }
  }
  npoints=j;

  /* end of dependences restoring */

  /* begin keep sources of same frequency */

  k=0;
  freq=Frequency(p);
  if (Nparents(p) > 0) {
    for (i=0;i<Nparents(p);i++) ppp[i]=(bodyptr) Parents(p)[i];
    for (i=0;i<Nparents(p);i++) {
      if (Frequency(ppp[i]) == freq) {
	Parents(p)[k]=(nodeptr) ppp[i];
	k++;
      }
    }
    Nparents(p)=k;
  }
	
  k=0;
  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;
    khelp=k;
  }

  k=0;
  freq=Frequency(p);
  if (Nsons(p) > 0) {
    for (i=0;i<Nsons(p);i++) sss[i]=(bodyptr) Sons(p)[i];
    for (i=0;i<Nsons(p);i++) {
      if (Frequency(sss[i]) == freq) {
	Sons(p)[k]=(nodeptr) sss[i];
	k++;
      }
    }
    Nsons(p)=k;
  }

  /* end keep sources of same frequency */

  check=0;

  npoints1=npoints;
  for (i=0;i<npoints;i++) { /* check for Flux or Fluxerror <= 0 */
    what1[i]=what[i];
    spectrum1[i]=spectrum[i];
  }
  
  k=0;
  for (i=0;i<npoints;i++) {
    if (Flux(spectrum1[i]) <= 0.0 || Fluxerror(spectrum1[i]) <= 0.0) spectrum1[i]=NULL;
    else k++;
  }
  npoints=k;
  for (i=0;i<npoints;i++) if (spectrum1[i] != NULL) {
    spectrum[i]=spectrum1[i];
    what[i]=what1[i];
  }

  npoints1=npoints;
  for (i=0;i<npoints;i++) { /* remember spectrum for second try with b=-0.9 */
    what1[i]=what[i];
    spectrum1[i]=spectrum[i];
    spectrum2[i]=spectrum[i];
  }


 jumpy:

  if (specrun > 0) {            /* restore remembered spectrum */
    npoints=npoints1;
    for (i=0;i<npoints;i++) {
      what[i]=what1[i];
      spectrum[i]=spectrum1[i];
    }
  }

  qq=sorting;                    /* sort spectrum according to frequency */
  for (i=0;i<npoints;i++) {
    Compnum(qq)=spectrum[i];
    Compfrequency(qq)=Frequency(spectrum[i]);
    Compi(qq)=what[i];
    qq++;
  }
  qsort(sorting,npoints,sizeof(comppare),ffcompare);

  check=check2=qwe=0;
  while(check == 0) {
    pl=1;
    qq=sorting;
    freqhelp[0]=Frequency(Compnum(qq));
    for (i=0;i<npoints;i++) {      /* count sources of different frequency */
      j=0;
      for (kk=0;kk<pl;kk++) if (Frequency(Compnum(qq)) == freqhelp[kk]) j=1;
      if (j == 0) {
	freqhelp[pl]=Frequency(Compnum(qq));
	pl++;
      }
      qq++;
    }

    if (pl > 2) {       /* if more than two sources at different frequencies */
      for (i=1;i<=npoints;i++) spflag[i]=0;
      qq=sorting;
      for (i=1;i<=npoints;i++) {
	nu[i]=Frequency(Compnum(qq));
	flux[i]=Flux(Compnum(qq));
	lognu[i]=log10(nu[i]);
	logflux[i]=log10(flux[i]);
	fluxerror[i]=Fluxerror(Compnum(qq));
	wat[i]=Compi(qq);
	qq++;
      }
      
      medfit(lognu,logflux,npoints,&a,&b,&abdev);  /* least absolute deviation fit */

    anfang:

      check1=0;
      for (i=1;i<=npoints;i++) { /* check for sources that fit into spectrum */
	helpp=a+b*lognu[i];      
	estflux[i]=pow(10.,helpp); 
	if (estflux[i] <= flux[i]+factor*fluxerror[i] &&  estflux[i] >= flux[i]-factor*fluxerror[i]) {
	  spflag[i]=1;
	  check1++;
	}
      }
      check4=0;

      /*printf("check1=%i  npoints=%i \n",check1,npoints);*/

      if (check1 == npoints) {
	check=npoints;
	check1=1;
	check4=1;
	goto ende;
      }
      if (check1 == npoints) check=1;
      if (check1 < check2) {
	a=aold;
	b=bold;
	zaehler++;
	if (zaehler == 1) goto anfang;
      } else {
	check2=check1;
	aold=a;
	bold=b;
      }
      zaehler=0;

      if (specrun > 0 && check1 > check2 && zaehler1 == 0) {
	specrun=0;
	zaehler1=1;
      }

      if (specrun > 0 && qwe < (npoints1-3)) {
	qwe++;
	for (i=1;i<=npoints;i++) spflag[i]=0;
	check1=0;
	for (i=1;i<=npoints;i++) { /* check for sources that fit into spectrum */
	  helpp=acalc+bcalc*lognu[i];      
	  estflux[i]=pow(10.,helpp); 
	  if (estflux[i] <= flux[i]+factor*fluxerror[i] &&  estflux[i] >= flux[i]-factor*fluxerror[i]) {
	    spflag[i]=1;
	    check1++;
	  }
	}
	if (check1 == npoints) check=1;
      }

      ps=0;
      helpq=0;
      for (i=1;i<=npoints;i++) {/* search for double frequency source with greatest deviation */
	helph=0;
	for (j=1;j<=npoints;j++) if (i != j && nu[i] == nu[j]) helph=1;
	helpp=fabs(estflux[i]-flux[i]);
	if (spflag[i] == 0 && wat[i] != SOURCE && helpp > helpq && helph == 1) {
	  helpq=helpp;
	  ps=i;
	}
      }
      if (ps == 0) {
	helpq=0;
	for (i=1;i<=npoints;i++) {/* search for source with greatest deviation */
	  helpp=fabs(estflux[i]-flux[i]);
	  if (spflag[i] == 0 && wat[i] != SOURCE && helpp > helpq) {
	    helpq=helpp;
	    ps=i;
	  }
	}
      }
      if (ps > 0) {
	spflag[ps]=-1;               /* flag source with greatest deviation */
	qq=sorting;
	for (i=0;i<npoints;i++) {
	  spectrum[i]=Compnum(qq);
	  what[i]=Compi(qq);
	  qq++;
	} 
	qq=sorting;
	k=0;
	for (i=0;i<npoints;i++) {   /* remove source with greatest deviation */
	  if (spflag[i+1] >= 0) {
	    Compnum(qq)=spectrum[i];
	    Compi(qq)=what[i];
	    qq++;
	    k++;
	  } 
	}
	if (npoints != k+1) {
	  printf("problem in spectra \n");
	  exit(0);
	}
	npoints=k;
      } else check=1;
    } else check=1;
  }              /* makes another try to fit a spectrum with one point less */

  check=0;
  pl=1;
  qq=sorting;
  freqhelp[0]=Frequency(Compnum(qq));
  for (i=0;i<npoints;i++) {      /* count sources of different frequency */
    j=0;
    for (kk=0;kk<pl;kk++) if (Frequency(Compnum(qq)) == freqhelp[kk]) j=1;
    if (j == 0) {
      freqhelp[pl]=Frequency(Compnum(qq));
      pl++;
    }
    qq++;
  }
  
  if (pl > 2 && check == 1) { /* if more than 2 independent point: final fit */
    qq=sorting;
    for (i=1;i<=npoints;i++) {
      nu[i]=Frequency(Compnum(qq));
      flux[i]=Flux(Compnum(qq));
      lognu[i]=log10(nu[i]);
      logflux[i]=log10(flux[i]);
      fluxerror[i]=Fluxerror(Compnum(qq));
      if (Wflag(Compnum(qq)) == 2) Resolved(Compnum(qq))=(int) Frequency(Compnum(qq));
      qq++;
    }
    medfit(lognu,logflux,npoints,&a,&b,&abdev);
    
    if (aold != a) {
      a=aold;
      b=bold;
    }
    check1=check=0;
    for (i=1;i<=npoints;i++) { /* check if initial source fits into spectrum */
      helpp=a+b*lognu[i];      
      estflux[i]=pow(10.,helpp);
      spflag[i]=0;
      if (estflux[i] <= flux[i]+factor*fluxerror[i] &&  estflux[i] >= flux[i]-factor*fluxerror[i]) {
	spflag[i]=1;
	check++;
	if (wat[i] == SOURCE) check1=1;
      }
    }

    pll=0;                        /* apply proximity criterion */
    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) check1=0;

    helpp=0.0;
    helpq=1.e10;
    for (i=0;i<npoints;i++) {
      if (Majaxis(spectrum[i]) == 0) helph=Resolution(spectrum[i]);
      else helph=Majaxis(spectrum[i]);
      helpp=MAX(helpp,helph);
      helpq=MIN(helpq,helph);
    }
    if (helpp/helpq > 10. && b < -2.0 && npoints == 3) check1=0;

  ende:
    check3=0;
    if (check == npoints && check1 == 1 && npoints >= 3 /*&& b <= -0.5*/) { /* if spectrum steeper than -0.5 */
      check3=1;
      nspectra++;
      k=Nbrothers(pskip);
      qq=sorting;
      for (i=1;i<=npoints;i++) {   /* keep sources at different frequencies as brothers */
	if (Compi(qq) != SOURCE && spflag[i] == 1) {
	  Brothers(pskip)[k]=(nodeptr) Compnum(qq);
	  k++;
	}
	qq++;
      }
      Nbrothers(pskip)=k;
      Specindx(p)=b;
      Absz(p)=a;
      if (specrun == 0) npoints2=npoints;
      if (check4 == 1) goto vorbei;
    }
  }

  specrun++;
  if (specrun == 1 && check3 == 0) {
    bcalc=-0.9;
    acalc=log10(Flux(p))-bcalc*log10(Frequency(p));
    goto jumpy;
  }
  
  if (specrun == 2 && npoints2 >= npoints) {
    specrun=-1;
    goto jumpy;
  }

  if (Specindx(p) == 0.0 && npoints1 > 2) {   /* take only maximum values if more than one point per frequency */
    for (i=0;i<npoints1;i++) if (Flux(spectrum1[i]) <= 0.0 || Fluxerror(spectrum1[i]) <= 0.0) spectrum1[i]=NULL;
    for (i=0;i<npoints1;i++) {
      for (j=0;j<npoints1;j++) {
	if (i != j && spectrum1[i] != NULL && spectrum1[j] != NULL) {
	  if (Frequency(spectrum1[i]) == Frequency(spectrum1[j])) {
	    if (Flux(spectrum1[i]) >= Flux(spectrum1[j])) spectrum1[j]=NULL;
	    else spectrum1[i]=NULL;
	  }
	}
      }
    }
    j=0;
    for (i=0;i<npoints1;i++) {
      if (spectrum1[j] != NULL) {
	spectrum[j]=spectrum1[i]; 
	j++;
      }
    }
    npoints=j;


    if (npoints > 2) {/* begin of least square fit */

      for (i=1;i<=npoints;i++) {
	nu[i]=Frequency(spectrum[i-1]);
	flux[i]=Flux(spectrum[i-1]);
	lognu[i]=log10(nu[i]);
	logflux[i]=log10(flux[i]);
	fluxerror[i]=Fluxerror(spectrum[i-1]);
      }
      
      k=1;
      fit(lognu,logflux,npoints,fluxerror,k,&a,&b,&siga,&sigb,&chi2,&sum);
      
      for (i=1;i<=npoints;i++) spflag[i]=0;
      check1=i=1;
      helpp=a+b*lognu[i];      
      estflux[i]=pow(10.,helpp);
      if (estflux[i] <= flux[i]+factor*fluxerror[i] &&  estflux[i] >= flux[i]-factor*fluxerror[i]) {
	spflag[i]=1;
	for (i=2;i<=npoints;i++) { /* check for sources that fit into spectrum */
	  helpp=a+b*lognu[i];      
	  estflux[i]=pow(10.,helpp); 
	  if (estflux[i] <= flux[i]+factor*fluxerror[i] &&  estflux[i] >= flux[i]-factor*fluxerror[i]) {
	    spflag[i]=1;
	    check1++;
	  }
	}

	pll=0;               /* apply proximity criterion */
	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) check1=0;

	helpp=0.0;
	helpq=1.e10;
	for (i=0;i<npoints;i++) {
	  if (Majaxis(spectrum[i]) == 0) helph=Resolution(spectrum[i]);
	  else helph=Majaxis(spectrum[i]);
	  helpp=MAX(helpp,helph);
	  helpq=MIN(helpq,helph);
	}
	if (helpp/helpq > 10. && b < -2.0 && npoints == 3) check1=0;

	if (check1 > 2) {
	  k=Nbrothers(pskip);
	  for (i=1;i<=npoints;i++) {   /* keep sources at different frequencies as brothers */
	    if (spectrum[i-1] != pskip && spflag[i] == 1) {
	      Brothers(pskip)[k]=(nodeptr) spectrum[i-1];
	      k++;
	    }
	  }
	  Nbrothers(p)=k;
	  Specindx(p)=b;
	  Absz(p)=a;
	} 
      }
    }
  }

  if (Specindx(p) == 0.0 && npoints1 > 2) {
    sprintf(Datei,"Abfall%i.dat",1);
    datei_ptr=fopen(Datei,"a");
    if (Pos(p)[0] >= em1 && Pos(p)[0] < em2) {
      j=1;
      for (i=0;i<npoints1;i++) if (Frequency(spectrum2[i]) != Frequency(p)) j++;
      if (j > 2 ) {
	fprintf(datei_ptr,"%s \n",Name(p));
	fprintf(datei_ptr,"%i  %10.8e %10.8e \n",j,Absz(p),Specindx(p));
	fprintf(datei_ptr,"%10.8e  %10.8e  %10.8e  \n",Frequency(p),Flux(p),Fluxerror(p));
	for (i=0;i<npoints1;i++) {
	  if (Frequency(spectrum2[i]) != Frequency(p)) {
	    pq=(bodyptr) spectrum2[i];
	    fprintf(datei_ptr,"%s\n",Name(pq));
	    fprintf(datei_ptr,"%10.8e   %10.8e   %10.8e  \n",Frequency(pq),Flux(pq),Fluxerror(pq));
	  }
	}
      }
    }
    fclose(datei_ptr);
  }

 vorbei:


  free_vector0(nu,1,maxspec);
  free_vector0(flux,1,maxspec);
  free_vector0(fluxerror,1,maxspec);
  free_vector0(logflux,1,maxspec);
  free_vector0(lognu,1,maxspec);
  free_vector0(estflux,1,maxspec);
  free_ivector0(spflag,1,maxspec);
  free_ivector0(wat,1,maxspec);
}



/*
 *  quicksort subroutine
 */ 
	
	 
local int ffcompare(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);
}