# include "dvo.h"

void cprecess (Average *average, int Naverage, double in_epoch, double out_epoch) {

  int i;
  double T;
  double A, D, RA, DEC, zeta, z, theta;
  double SA, CA, SD, CD;
  
  T = (out_epoch - in_epoch) / 100.0;
  
  zeta  = RAD_DEG*(0.6406161*T + 0.0000839*T*T + 0.0000050*T*T*T);
  theta = RAD_DEG*(0.5567530*T - 0.0001185*T*T - 0.0000116*T*T*T);
  z     =          0.6406161*T + 0.0003041*T*T + 0.0000051*T*T*T;
  
  for (i = 0; i < Naverage; i++) {
    A = average[i].R;
    D = average[i].D;
    SD =  cos(RAD_DEG*A + zeta)*sin(theta)*cos(RAD_DEG*D) + cos(theta)*sin(RAD_DEG*D);
    CD = sqrt (1 - SD*SD);
    SA =  sin(RAD_DEG*A + zeta)*cos(RAD_DEG*D)/CD;
    CA = (cos(RAD_DEG*A + zeta)*cos(theta)*cos(RAD_DEG*D) - sin(theta)*sin(RAD_DEG*D))/CD;
    
    DEC = DEG_RAD*asin(SD);
    RA  = DEG_RAD*atan2(SA, CA) + z;
    
    if (RA < 0)
      RA += 360;
    
    average[i].R = RA;
    average[i].D = DEC; 
  }

}

/* values are ave[i].R, ave[i].D, ave[i].M */
void sortave (Average *ave, int N) {

  int l,j,ir,i;
  Average tmp;

  if (N < 2) return;
  l = N >> 1;
  ir = N - 1;
  for (;;) {
    if (l > 0) {
      l--;
      tmp = ave[l];
    }
    else {
      tmp = ave[ir];
      ave[ir] = ave[0];
      if (--ir == 0) {
	ave[0] = tmp;
	return;
      }
    }
    i = l;
    j = (l << 1) + 1;
    while (j <= ir) {
      if (j < ir && ave[j].R < ave[j+1].R) j++;
      if (tmp.R < ave[j].R) {
	ave[i] = ave[j];
	j += (i=j) + 1;
      }
      else j = ir + 1;
    }
    ave[i] = tmp;
  }
}

/* values are ave[i].R, ave[i].D, ave[i].M */
void sort_images (Image *image, int N) {

  int l,j,ir,i;
  Image tmp;

  if (N < 2) return;
  l = N >> 1;
  ir = N - 1;
  for (;;) {
    if (l > 0) {
      l--;
      tmp = image[l];
    }
    else {
      tmp = image[ir];
      image[ir] = image[0];
      if (--ir == 0) {
	image[0] = tmp;
	return;
      }
    }
    i = l;
    j = (l << 1) + 1;
    while (j <= ir) {
      if (j < ir && image[j].tzero < image[j+1].tzero) j++;
      if (tmp.tzero < image[j].tzero) {
	image[i] = image[j];
	j += (i=j) + 1;
      }
      else j = ir + 1;
    }
    image[i] = tmp;
  }
}

/* sort subset by image[subset[i]].tzero */
void sort_image_subset (Image *image, int *subset, int N) {

  int l, j, ir, i, tmp;

  if (N < 2) return;
  l = N >> 1;
  ir = N - 1;
  for (;;) {
    if (l > 0) {
      l--;
      tmp = subset[l];
    }
    else {
      tmp = subset[ir];
      subset[ir] = subset[0];
      if (--ir == 0) {
	subset[0] = tmp;
	return;
      }
    }
    i = l;
    j = (l << 1) + 1;
    while (j <= ir) {
      if (j < ir && image[subset[j]].tzero < image[subset[j+1]].tzero) j++;
      if (image[tmp].tzero < image[subset[j]].tzero) {
	subset[i] = subset[j];
	j += (i=j) + 1;
      }
      else j = ir + 1;
    }
    subset[i] = tmp;
  }
}
