/*      FILE:  search.c
 *   PURPOSE: 
 *    Find spatial offset between the two images using a simple 
 *    cross-correlation technique.
 *    AUTHOR:  Kenneth J. Mighell (mighell@noao.edu)
 *  LANGUAGE:  ANSI C
 *      DATE:  2001JUL03
 * COPYRIGHT:  (C) 2001 Assoc. of Universities for Research in Astronomy Inc.
 */

#include "mx.h"
#include "qdphot.h"

#define DEBUG (MX_FALSE)

int qdphot5_Search_f3 (
  qdphot5_ParS *ParS,
  qdphot5_ImageS *ImageS1,
  qdphot5_ImageS *ImageS2
){

  char mxfunc[] = "qdphot5_Search_f3";
  int mxstatus;

  int status;
  int mxnstr;

  int nstr;
  int id;
  qdphot5_PhotS PhotS1;
  qdphot5_PhotS PhotS2;
  FILE *out1fd;
  FILE *out2fd;
  double x1cood;
  double y1cood;
  int rejects;
  int ok;
  int k;
  int j;
  int xi;
  int yi;
  double old_dxd;
  double old_dyd;
  int show;

#define XCORSIZE (10001)
  double bxs[XCORSIZE];
  double bys[XCORSIZE];
  double rxs[XCORSIZE];
  double rys[XCORSIZE];

  double big;
  double xcor;
  double dx;
  double dy;
  int radius;
  int jb;
  int je;
  int kb;
  int ke;

  mxnstr = ParS->mxstars;
  rejects = ParS->rejects;
  show = ParS->showsearch;
  out1fd = (FILE *)NULL;
  out2fd = (FILE *)NULL;

  /*
   * Initialize slices
   */
  for ( k=0; k<XCORSIZE; ++k ) {
    bxs[k] = 0.0;
    bys[k] = 0.0;
    rxs[k] = 0.0;
    rys[k] = 0.0;
  }

  if (show) printf( "\nSEARCH: begin ===================================\n\n");

/* 
 * 2nd image =================================================================
 */

  if (show) printf( "*** 2nd image ***\n" );

  if (show) printf( "%s\n", ImageS2->ifn );

  if (DEBUG) printf( "ParS->coodim is %d\n", ParS->coodim );

  ImageS2->dxd = 0.0;
  ImageS2->dyd = 0.0;
  status = qdphot5_ParS_CooS_Peaker_f3( ParS, ImageS2, 0 );
  if (status) {mxstatus = 1; goto mx_error1;}

  if (show) printf( " --> found %d objects\n", ParS->coostars );
  
  if (ParS->coostars<1) goto bye;

  for ( k=0; k<(ParS->coostars); ++k ) {

    /* coordinate file true X and Y */ 
    x1cood = ParS->CooS[k].xd;
    y1cood = ParS->CooS[k].yd;
    id = ParS->CooS[k].idi;

    qdphot5_PhotS_Calc_f5 (ImageS2, &PhotS2, x1cood, y1cood, id);

    if (PhotS2.ok) {
      xi = PhotS2.x1outd + 0.4999;
      if ( (xi>0) && (xi<XCORSIZE) ) rxs[xi] = rxs[xi] + PhotS2.sgnl;
      yi = PhotS2.y1outd + 0.4999;
      if ( (yi>0) && (yi<XCORSIZE) ) rys[yi] = rys[yi] + PhotS2.sgnl;
    }
      
  }

/* 
 * 1st image =================================================================
 */

  if (show) printf( "*** 1st image ***\n" );

  if (show) printf( "%s\n", ImageS1->ifn );

  if (DEBUG) printf( "ParS->coodim is %d\n", ParS->coodim );

  ImageS1->dxd = 0.0;
  ImageS1->dyd = 0.0;
  status = qdphot5_ParS_CooS_Peaker_f3( ParS, ImageS1, 0 );
  if (status) {mxstatus = 1; goto mx_error1;}

  if (show) printf( " --> found %d objects\n", ParS->coostars );
  
  if (ParS->coostars<1) goto bye;

  for ( k=0; k<(ParS->coostars); ++k ) {

    /* coordinate file true X and Y */ 
    x1cood = ParS->CooS[k].xd;
    y1cood = ParS->CooS[k].yd;
    id = ParS->CooS[k].idi;

    qdphot5_PhotS_Calc_f5 (ImageS1, &PhotS1, x1cood, y1cood, id);

    if (PhotS1.ok) {
      xi = PhotS1.x1outd + 0.4999;
      if ( (xi>0) && (xi<XCORSIZE) ) bxs[xi] = bxs[xi] + PhotS1.sgnl;
      yi = PhotS1.y1outd + 0.4999;
      if ( (yi>0) && (yi<XCORSIZE) ) bys[yi] = bys[yi] + PhotS1.sgnl;
    }
      
  }

  if (DEBUG) {
    printf( "\n" );
    for ( k=100; k<700; ++k ) {
      printf( "%d %10g %10g %10g %10g\n", k, bxs[k], bys[k], rxs[k], rys[k] );
    }
    printf( "\n" );
  }

/* 
 * Find shift in X ===========================================================
 */

  radius = ParS->search2 + 0.9999;
  kb = -radius;
  ke = radius;
  jb = radius + 1;
  je = (XCORSIZE-1) - (radius+1) ; /* 700 */

  if (show) printf( "\n Find shift in X\n" );
  big = 0.0;
  dx = -QDPHOT_BIG;
  /*  for ( k=-40; k<=40; ++k ) { */
  for ( k=kb; k<=ke; ++k ) {
    xcor = 0.0;
    /* for ( j=100; j<=700; ++j ) { */
    for ( j=jb; j<=je; ++j ) {
      xcor = xcor + bxs[j]*rxs[j+k];
    }
    if (xcor>big) {
      big = xcor;
      dx = k;
      if (show) printf( "%6d %14.1f <--\n", k, xcor );
    } else {
      if (show) printf( "%6d %14.1f    \n", k, xcor );
    }
  }
  if (show) printf( "\n DeltaX is %g\n", dx );

/* 
 * Find shift in Y ===========================================================
 */

  if (show) printf( "\n Find shift in Y\n" );
  big = 0.0; 
  dy = -QDPHOT_BIG;
  /*  for ( k=-40; k<=40; ++k ) { */
  for ( k=kb; k<=ke; ++k ) {
    xcor = 0.0;
    /* for ( j=100; j<=700; ++j ) { */
    for ( j=jb; j<=je; ++j ) {
      xcor = xcor + bys[j]*rys[j+k];
    }
    if (xcor>big) {
      big = xcor;
      dy = k;
      if (show) printf( "%6d %14.1f <--\n", k, xcor );
    } else {
      if (show) printf( "%6d %14.1f    \n", k, xcor );
    }
  }
  if (show) printf( "\n DeltaY is %g\n", dy );

  if (show) printf( "\n DeltaX is %8.2f  and DeltaY is %8.2f\n", dx, dy );

ok:
  mxstatus = 0;
  ParS->odx2 = dx;
  ParS->ody2 = dy;
  ImageS2->dxd = dx;
  ImageS2->dyd = dy;
  ParS->dx[1] = ParS->dx[0];
  ParS->dx_is_INDEF[1] = MX_TRUE;
  ParS->dx[2] = ParS->dx[0];
  ParS->dx_is_INDEF[2] = MX_TRUE;
  ParS->dy[1] = ParS->dy[0];
  ParS->dy_is_INDEF[1] = MX_TRUE;
  ParS->dy[2] = ParS->dy[0];
  ParS->dy_is_INDEF[2] = MX_TRUE;

  if (show) printf( "\nSEARCH: end =====================================\n\n");

  goto bye;

mx_error1:
  mxp_tmpmsg_init_f0();
mx_error2:
  mxp_errmsg_append_f3 (mxfunc, mxstatus, MX.tmpmsg);
  goto bye;

bye:
  return(mxstatus);

}
#undef DEBUG
/* end-of-file */