/*

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                   !!
!!                   GNU General Public License                      !!
!!                                                                   !!
!! This file is part of the Flexible Modeling System (FMS).          !!
!!                                                                   !!
!! FMS is free software; you can redistribute it and/or modify       !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation.             !!
!!                                                                   !!
!! FMS is distributed in the hope that it will be useful,            !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of    !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the     !!
!! GNU General Public License for more details.                      !!
!!                                                                   !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to:                                 !!
!!          Free Software Foundation, Inc.                           !!
!!          59 Temple Place, Suite 330                               !!
!!          Boston, MA  02111-1307  USA                              !!
!! or see:                                                           !!
!!          http://www.gnu.org/licenses/gpl.txt                      !!
!!                                                                   !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

*/
#include <stdlib.h>
#include <math.h>
#include "mosaic_util.h"
#include "get_contact.h"


/************************************************************************************************
         int get_align_contact(
  This routine will return number of algined contacts bewteen two tiles (line contact).
  This routine assume the starting and ending points of the  contact line are coincidence with
  the grid points of both tiles. lrg_rectangle tiles are assumed. 
  It will return the contact information, which includes 

*************************************************************************************************/

double* east_bound(const double *data, int nx, int ny);
double* west_bound(const double *data, int nx, int ny);
double* south_bound(const double *data, int nx, int ny);
double* north_bound(const double *data, int nx, int ny);
#define EPSLN (1.0e-10)

int get_contact_index( int size1, int size2, double *x1, double *y1, double *x2, double *y2, double periodx,
		       double periody, int *start1, int *end1, int *start2, int *end2);

int get_align_contact(int tile1, int tile2, int nx1, int ny1, int nx2, int ny2, 
                      const double *x1, const double *y1, const double *x2, 
                      const double *y2, double periodx, double periody,
		      int *istart1, int *iend1, int *jstart1, int *jend1, 
                      int *istart2, int *iend2, int *jstart2, int *jend2)
{
  double *xb1, *yb1, *xb2, *yb2;
  int ncontact, start1, end1, start2, end2;
  
  ncontact = 0;
  /* East bound of tile1 and west bound of tile2 */
  xb1 = east_bound(x1, nx1, ny1);
  yb1 = east_bound(y1, nx1, ny1);
  xb2 = west_bound(x2, nx2, ny2);
  yb2 = west_bound(y2, nx2, ny2);
  if( get_contact_index( ny1, ny2, xb1, yb1, xb2, yb2, periodx, 0.0, &start1, &end1, &start2, &end2) ) {
    istart1[ncontact] = nx1-1;
    iend1[ncontact]   = nx1-1;
    istart2[ncontact] = 1;
    iend2[ncontact]   = 1;    
    jstart1[ncontact] = start1;
    jend1[ncontact]   = end1;
    jstart2[ncontact] = start2;
    jend2[ncontact]   = end2;
    ncontact++;   
  }

  /* East bound of tile1 and SOUTH bound of tile2, tile1 and tile must be different tile */
  if(tile1 != tile2) {
    free(xb2);
    free(yb2);
    xb2 = south_bound(x2, nx2, ny2);
    yb2 = south_bound(y2, nx2, ny2);
    if( get_contact_index( ny1, nx2, xb1, yb1, xb2, yb2, 0.0, 0.0, &start1, &end1, &start2, &end2) ) {
      istart1[ncontact] = nx1-1;
      iend1[ncontact]   = nx1-1;
      istart2[ncontact] = start2;
      iend2[ncontact]   = end2;    
      jstart1[ncontact] = start1;
      jend1[ncontact]   = end1;
      jstart2[ncontact] = 1;
      jend2[ncontact]   = 1;
      ncontact++;
    }
  }
  free(xb1);
  free(yb1);   
  free(xb2);
  free(yb2);  
  
  /* South bound of tile1 and NORTH bound of tile2 */
  xb1 = south_bound(x1, nx1, ny1);
  yb1 = south_bound(y1, nx1, ny1);
  xb2 = north_bound(x2, nx2, ny2);
  yb2 = north_bound(y2, nx2, ny2);
  if( get_contact_index( nx1, nx2, xb1, yb1, xb2, yb2, 0.0, periody, &start1, &end1, &start2, &end2) ) {
    istart1[ncontact] = start1;
    iend1[ncontact]   = end1;
    istart2[ncontact] = start2;
    iend2[ncontact]   = end2;    
    jstart1[ncontact] = 1;
    jend1[ncontact]   = 1;
    jstart2[ncontact] = ny2-1;
    jend2[ncontact]   = ny2-1;
    ncontact++;
  }

  /* South bound of tile1 and East bound of tile2, tile1 and tile must be different tile*/  
  if(tile1 != tile2 ) {
    free(xb2);
    free(yb2);
    xb2 = east_bound(x2, nx2, ny2);
    yb2 = east_bound(y2, nx2, ny2);
    if( get_contact_index( nx1, ny2, xb1, yb1, xb2, yb2, 0.0, 0.0, &start1, &end1, &start2, &end2) ) {
      istart1[ncontact] = start1;
      iend1[ncontact]   = end1;
      istart2[ncontact] = nx2-1;
      iend2[ncontact]   = nx2-1;    
      jstart1[ncontact] = 1;
      jend1[ncontact]   = 1;
      jstart2[ncontact] = start2;
      jend2[ncontact]   = end2;
      ncontact++;
    }
  }
  free(xb1);
  free(yb1);   
  free(xb2);
  free(yb2);  

  /* to avoid duplicate, the following will be done only when tile1 not equal to tile2 */
  if(tile1 != tile2) {
    /* West bound of tile1 and east bound of tile2*/
    xb1 = west_bound(x1, nx1, ny1);
    yb1 = west_bound(y1, nx1, ny1);
    xb2 = east_bound(x2, nx2, ny2);
    yb2 = east_bound(y2, nx2, ny2);
    if( get_contact_index( ny1, ny2, xb1, yb1, xb2, yb2, periodx, 0.0, &start1, &end1, &start2, &end2) ) {
      istart1[ncontact] = 1;
      iend1[ncontact]   = 1;
      istart2[ncontact] = nx2-1;
      iend2[ncontact]   = nx2-1;    
      jstart1[ncontact] = start1;
      jend1[ncontact]   = end1;
      jstart2[ncontact] = start2;
      jend2[ncontact]   = end2;
      ncontact++;
    }
    free(xb2);
    free(yb2);

    /* West bound of tile1 and North bound of tile2 */  
    xb2 = north_bound(x2, nx2, ny2);
    yb2 = north_bound(y2, nx2, ny2);
    if( get_contact_index( ny1, nx2, xb1, yb1, xb2, yb2, 0.0, 0.0, &start1, &end1, &start2, &end2) ) {
      istart1[ncontact] = 1;
      iend1[ncontact]   = 1;
      istart2[ncontact] = start2;
      iend2[ncontact]   = end2;    
      jstart1[ncontact] = start1;
      jend1[ncontact]   = end1;
      jstart2[ncontact] = ny2-1;
      jend2[ncontact]   = ny2-1;
      ncontact++;
    }
    free(xb1);
    free(yb1);   
    free(xb2);
    free(yb2);  

  
    /* North bound of tile1 and South bound of tile2 */
    xb1 = north_bound(x1, nx1, ny1);
    yb1 = north_bound(y1, nx1, ny1);
    xb2 = south_bound(x2, nx2, ny2);
    yb2 = south_bound(y2, nx2, ny2);
    if( get_contact_index( nx1, nx2, xb1, yb1, xb2, yb2, 0.0, periody, &start1, &end1, &start2, &end2) ) {
      istart1[ncontact] = start1;
      iend1[ncontact]   = end1;
      istart2[ncontact] = start2;
      iend2[ncontact]   = end2;    
      jstart1[ncontact] = ny1-1;
      jend1[ncontact]   = ny1-1;
      jstart2[ncontact] = 1;
      jend2[ncontact]   = 1;
      ncontact++;
    }
    free(xb2);
    free(yb2);

    /* North bound of tile1 and West bound of tile2 */  
    xb2 = west_bound(x2, nx2, ny2);
    yb2 = west_bound(y2, nx2, ny2);
    if( get_contact_index( nx1, ny2, xb1, yb1, xb2, yb2, 0.0, 0.0, &start1, &end1, &start2, &end2) ) {
      istart1[ncontact] = start1;
      iend1[ncontact]   = end1;
      istart2[ncontact] = 1;
      iend2[ncontact]   = 1;    
      jstart1[ncontact] = ny1-1;
      jend1[ncontact]   = ny1-1;
      jstart2[ncontact] = start2;
      jend2[ncontact]   = end2;
      ncontact++;
    }
    free(xb1);
    free(yb1);   
    free(xb2);
    free(yb2);  
  }

  /* when tile1 = tile2, we need to consider about folded. Only foled north is considered here */
  if(tile1 == tile2) {
    int i, folded = 1;
    double dx;
    xb1 = north_bound(x1, nx1, ny1);
    yb1 = north_bound(y1, nx1, ny1);
    for(i=0; i<nx1/2; i++) {
      if( yb1[i] != yb1[nx1-i-1] ) {
	folded = 0;
	break;
      }
      dx = fabs(xb1[i] - xb1[nx1-i-1]);
      if( dx !=0 && dx != 360 ) {
	folded = 0;
	break;
      }
    }
    if(folded) {
      istart1[ncontact] = 1;
      iend1[ncontact]   = nx1/2;
      istart2[ncontact] = nx1-1;
      iend2[ncontact]   = nx1/2+1;    
      jstart1[ncontact] = ny1-1;
      jend1[ncontact]   = ny1-1;
      jstart2[ncontact] = ny1-1;
      jend2[ncontact]   = ny1-1;
      ncontact++;
    }
    free(xb1);
    free(yb1);
  }
  
  return ncontact;
  
};


int get_contact_index( int size1, int size2, double *x1, double *y1, double *x2, double *y2, double periodx,
		       double periody, int *start1, int *end1, int *start2, int *end2)
{
  int i1, i2;
  double dx, dy;

  /* Find the first point in tile 1 cocindent with a point in tile2  */
  *start1 = -1;
  *start2 = -1;
  for(i1=0; i1<size1; i1++) {
    for(i2=0; i2<size2; i2++) {
      dx = fabs(x1[i1]- x2[i2]);
      dx = min(dx, fabs(dx-periodx));
      dy = fabs(y1[i1]- y2[i2]);
      dy = min(dy, fabs(dy-periody));
      if( dx < EPSLN && dy <EPSLN ) {
	*start1 = i1+1;
	*start2 = i2+1;
	goto foundstart;
      }
    }
  }

  return 0;
    
  foundstart:

  /* Find the last point in tile 1 cocindent with a point in tile2 */
  *end1 = -1;
  *end2 = -1;
  for(i1=size1-1; i1>=0; i1--) {
    for(i2=size2-1; i2>=0; i2--) {
      dx = fabs(x1[i1]- x2[i2]);
      dx = min(dx, fabs(dx-periodx));
      dy = fabs(y1[i1]- y2[i2]);
      dy = min(dy, fabs(dy-periody));
      if( dx < EPSLN && dy <EPSLN ) {
	*end1 = i1+1;
	*end2 = i2+1;
	goto foundend;
      }
    }
  }

  return 0;
    
 foundend: if( *start1 == *end1 || *start2 == *end2 ) return 0;

  if(*start1 > *end1 )
    (*start1)--;
  else
    (*end1)--;

  if(*start2 > *end2 )
    (*start2)--;
  else
    (*end2)--;  

  return 1;
    
};


double* west_bound(const double *data, int nx, int ny)
{
  int i;
  double *bnd;
  bnd = (double *)malloc(ny*sizeof(double));
  for(i=0; i<ny; i++) bnd[i] = data[i*nx];
  return bnd;
}

double* east_bound(const double *data, int nx, int ny)
{
  int i;
  double *bnd;
  bnd = (double *)malloc(ny*sizeof(double));
  for(i=0; i<ny; i++) bnd[i] = data[i*nx+nx-1];
  return bnd;
}

double* south_bound(const double *data, int nx, int ny)
{
  int i;
  double *bnd;
  bnd = (double *)malloc(nx*sizeof(double));
  for(i=0; i<nx; i++) bnd[i] = data[i];
  return bnd;
}

double* north_bound(const double *data, int nx, int ny)
{
  int i;
  double *bnd;
  bnd = (double *)malloc(nx*sizeof(double));
  for(i=0; i<nx; i++) bnd[i] = data[(ny-1)*nx+i];
  return bnd;
}