#include <string.h>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>

/**
 * The replace function
 *
 * Searches all of the occurrences using recursion
 * and replaces with the given string
 * @param char * o_string The original string
 * @param char * s_string The string to search for
 * @param char * r_string The replace string
 * @return void The o_string passed is modified
 */
void replace(char * o_string, char * s_string, char * r_string) {
      //a buffer variable to do all replace things
      char buffer[255];
      //to store the pointer returned from strstr
      char * ch;

      //first exit condition
      if(!(ch = strstr(o_string, s_string)))
              return;

      //copy all the content to buffer before the first occurrence of the search string
      strncpy(buffer, o_string, ch-o_string);

      //prepare the buffer for appending by adding a null to the end of it
      buffer[ch-o_string] = 0;

      //append using sprintf function
      sprintf(buffer+(ch - o_string), "%s%s", r_string, ch + strlen(s_string));

      //empty o_string for copying
      o_string[0] = 0;
      strcpy(o_string, buffer);
      //pass recursively to replace other occurrences
      return replace(o_string, s_string, r_string);
 }



void SortedUniqueElements(float in_array[], float out_array[], int size, int *NU)
{
   int i, j;
   float temp;
   *NU=0;
   //printf("size: %d\n",size);
   for (i = 0; i < size; i++)
     {
       for (j = i+1; j < size; j++) if (in_array[i] == in_array[j]) break;
       if (j == size) out_array[(*NU)++] = in_array[i];
     }
   //printf("NU: %d\n",*NU);

   for ( i = 0; i < (*NU)-1; i++)
     {
       for ( j = i+1; j < *NU; j++)
	 {
	   if (out_array[i] > out_array[j])
	     {
	       temp = out_array[i];
	       out_array[i] = out_array[j];
	       out_array[j] = temp;
	     }
	 }
     }
}



int  interp1D (float x, float *v, int *n, int nmax)
{
  int j, iok; //float tmp;
  for (j=0, iok=0; j<*n; j++) 
    if (x==*(v+j)) // if value already exists
      { iok=1; break; }
  if (iok==0) // it does not exist, appends value
    { *(v+j) = x ; (*n)++; }
  if (*n>nmax)
    {
      printf("Error: maximum in dusty table (%d/%d) reached\n", *n, nmax);
      for (j=0; j<nmax; j++) {printf(" %.3g",*(v+j));} printf("\n");
      exit(1);
    }
  return j; //returns assigned index
}



// return interval (i1,i2) and interpolation factor fact on a 1D vector *v
// aditional parameters: 
// flagrange= -1 for only negative values, +1 for only positive, 0 for all
//    (but it requires negative/positive intervals to be contiguous)
// flaglimits= 0 for no-extrapolation (stay within limits i1,i2), 
//     -1 for extrapolation in the lower limit
//     1 for extrapolation in the upper limit
//     2 for extrapolation in both limits
float locateitin1Dtab (float x, float *v, int n, int *i1, int *i2,
		       int flagrange, int flaglimits, int flaglog)
{
  int j, lim1, lim2, nlim; 
  float fact;
  //linear: flaglog=0; logarithmic: flaglog=1
  //absolute limits: (lim1,lim2,nlim) instead of (0,n-1,n)
  lim1=0; lim2=n-1; nlim=n;
  for (j=0; j<n; j++)
    {
      if (flagrange==-1 && *(v+j)>0.0) continue;
      if (flagrange==+1 && *(v+j)<0.0) continue;
      lim1=j; break; // first point in interval
    }
  for (j=n-1; j>=0; j--)
    {
      if (flagrange==-1 && *(v+j)>0.0) continue;
      if (flagrange==+1 && *(v+j)<0.0) continue;
      lim2=j; break; // last point in interval
    }
  nlim = lim2-lim1;
  // now search within these limits:
  if (nlim==1)  // || x <= *(v+lim1)) //PROBLEM HERE?
    {*i1=lim1; *i2=lim2;}
  else if  (x == *(v+lim2)) //giada 30/01/2018
    {*i1=lim2-1; *i2=lim2; }
  else if (x < *(v+lim1)) //beyond lower limit 
    {*i1=lim1; *i2=lim1+1;}
  else if (x > *(v+lim2)) //beyond upper limit 
    {*i1=lim2-1; *i2=lim2;}
  else
    for (j=lim1; j<lim2; j++)
      {
	if ( (x >= *(v+j) && x < *(v+j+1)) || //this accepts non-monotonic tables:
	     (x <= *(v+j) && x > *(v+j+1)) )
	{ *i1=j; *i2=j+1; break; }
     }

  if (flaglog == 1)
     fact = (*i1==*i2) ? 1.0 :
       (log10(x) - log10(*(v+*i2))) / (log10(*(v+*i1)) - log10(*(v+*i2)));
  else if (flaglog == 0)
     fact = (*i1==*i2) ? 1.0 :
       (x - *(v+*i2)) / (*(v+*i1) - *(v+*i2));
  else{
     printf("wrong interpolation flag!\n");
     exit(1);}

  if (flaglimits==0) // avoids extrapolation
    {
      if (fact<0.0) fact=0.0;
      if (fact>1.0) fact=1.0;
    }
  else if (flaglimits==-1) // allows extrapolation to smaller values only 
    {if (fact>1.0) fact=1.0;}
  else if (flaglimits==+1) // allows extrapolation to larger values only
    if (fact<0.0) fact=0.0;
  // if flaglimits==2, nothing changes -> all extrapolations are allowed
  //printf ("%f %d %d %f\n", x, *i1, *i2, fact);
  return fact;
}