#include<interpola.h>
#include<stdlib.h>
#include<stdio.h>
#include<iostream.h>

void interpolad(int n,int m,  double *Y, double *X, int nr, double *x, double *res){
  int pos;
  int iu, im, il;
  bool ascnd;
  double coeff1, coeff2;
  ascnd = ( X[n-1] >= X[0] );
  
  for (int i=0; i<nr; i++) {
    
    il = -1;
    iu = n;
    while( iu-il > 1){
      im = (iu+il) >> 1;
      if ( x[i]  >= X[im] == ascnd )
	il = im;
      else
	iu = im;
    }
    
    if( x[i]==X[0] ) pos = 0;
    else if ( x[i ]==X[n-1] ) pos = n-2;
    else pos = il;

    if( pos==-1) pos=0;
    if(pos==n-1) pos=n-2;

    coeff1 = ( x[i]    -  X[pos]  )/(X[pos+1]-X[pos]);
    coeff2 = ( X[pos+1]  -  x[i]  )/(X[pos+1]-X[pos]);
    for(int j=0; j<m; j++) {
      res[i*m +j] = coeff1*Y[ j + (pos+1)*m ]  + coeff2*Y[  j+ pos*m ];
    }
  }
}

void interpolaD(int n,int m,  complex<double> *Y, double *X, int nr, double *x, complex<double> *res){
  int pos;
  int iu, im, il;
  bool ascnd;
  double coeff1, coeff2;
  ascnd = ( X[n-1] >= X[0] );
  
  for (int i=0; i<nr; i++) {
    
    il = -1;
    iu = n;
    while( iu-il > 1){
      im = (iu+il) >> 1;
      if ( x[nr]  >= X[im] == ascnd )
	il = im;
      else
	iu = im;
    }
    
    if( x[nr]==X[0] ) pos = 0;
    else if ( x[nr]==X[n-1] ) pos = n-2;
    else pos = il;
    coeff1 = ( x[nr]    -  X[pos]  )/(X[pos+1]-X[pos]);
    coeff2 = ( X[pos+1]  -  x[nr]  )/(X[pos+1]-X[pos]);
    for(int j=0; j<m; j++) {
      res[i*m +j] = coeff1*Y[ i*m + pos+1 ]  + coeff2*Y[  i*m + pos ];
    }
  }
}
