/* Last edited on 2008-01-22 by ptpinho */ 
/* This file contains the functions to compute derivatives*/

 /* We need to define _GNU_SOURCE because we use the constant NAN. */
#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>


#include "Int2Dlocal.h"
#include "deriv.h"


 /*==============================================================================================*/
    /* Cálculo da derivada de uma função segundo x ou y (flag=1 ou flag=0) num ponto ix ou iy */
                                 /* DERIVADA ANTI-SIMETRICA */
 /*==============================================================================================*/
 double DerivAS(double **u, double dxy, int ix, int jy, int N, int flag)
 { 
   int P = 4;           /* Number of interpolation points. */
   int sk[P];	        /* Indices of the samples needed to compute the derivative */
   double h[P];	        /* Weights needed to compute the derivative */
   double v[P];	        /* Function values needed to compute the derivative */
   int i, k;
   double du;

   if (flag==1) k = ix; else k = jy;    /*If  flag=1 then x derivative else y derivative */
	
   if (k==0)
     {
       sk[0] = 2;
       sk[1] = 1;
       sk[2] = 1;
       sk[3] = 2;
       h[0] = -1./12.;
       h[1] = 2./3.;
       h[2] = 2./3.;
       h[3] = -1./12.;
     }
   else if (k==1)
     {
       sk[0] = 1;
       sk[1] = 0;
       sk[2] = 2;
       sk[3] = 3;
       h[0] = -1./12.;
       h[1] =-2./3.;
       h[2] = 2./3.;
       h[3] = -1./12.;
     }
   else if (k==N-1)
     {
       sk[0] = N-3;
       sk[1] = N-2;
       sk[2] = N;
       sk[3] = N-1;
       h[0] = 1./12.;
       h[1] =-2./3.;
       h[2] = 2./3.;
       h[3] = 1./12.;
     }
   else if (k==N)
     {
       sk[0] = N-2;
       sk[1] = N-1;
       sk[2] = N-1;
       sk[3] = N-2;
       h[0] = 1./12.;
       h[1] =-2./3.;
       h[2] =-2./3.;
       h[3] = 1./12.;
     }
   else
     {
       sk[0] = k-2;
       sk[1] = k-1;
       sk[2] = k+1;
       sk[3] = k+2;
       h[0] = 1./12.;
       h[1] =-2./3.;
       h[2] = 2./3.;
       h[3] =-1./12.;
     }
    /* compute the values of the function */
   for (i=0; i <P; i++)
     {
      /* if (sk[i] < 0 || sk[i] > N)
         {
           fprintf(stderr,"ERROR: deriv(): sk[%d]=%d\n", i, sk[i]);
           printf(" ERRO no calculo da derivada Anti-Simétrica\n");
           exit(-1);
         }*/
       if (flag==1)
	 v[i] = u[sk[i]][jy]; 
       else
	 v[i] = u[ix][sk[i]];     
     }
   /* compute the derivative */
   du = 0;
   for (i=0; i < P; i++) 
     du += v[i] * h[i];

   du /= (1*dxy);

   return du;
 }

 /*==============================================================================================*/
    /* Cálculo da derivada de uma função segundo x ou y (flag=1 ou flag=0) num ponto ix ou iy */
                                      /* DERIVADA SIMETRICA */
 /*==============================================================================================*/

double DerivS(double **u, double dxy, int ix, int jy, int N, int flag)
 { 
   int P = 4;         /* Number of interpolation points. */
   int sk[P];	      /* the positions of the samples needed to compute the derivative */
   double h[P];	      /* the weights needed to compute the derivative */
   double v[P];	      /* the function values needed to compute the derivative */
   int i, k;
   double du;

   if (flag==1) k = ix; else k = jy;   /*If  flag=1 then x derivative else y derivative */
  
   if (k==0)
     {
       sk[0] = 2;
       sk[1] = 1;
       sk[2] = 1;
       sk[3] = 2;
       h[0] = 1./12.;
       h[1] = -2./3.;
       h[2] = 2./3.;
       h[3] = -1./12.;
     }
   else if (k==1)
     {
       sk[0] = 1;
       sk[1] = 0;
       sk[2] = 2;
       sk[3] = 3;
       h[0] = 1./12.;
       h[1] =-2./3.;
       h[2] = 2./3.;
       h[3] = -1./12.;
     }
   else if (k==N-1)
     {
       sk[0] = N-3;
       sk[1] = N-2;
       sk[2] = N;
       sk[3] = N-1;
       h[0] = 1./12.;
       h[1] =-2./3.;
       h[2] = 2./3.;
       h[3] = -1./12.;
     }
   else if (k==N)
     {
       sk[0] = N-2;
       sk[1] = N-1;
       sk[2] = N-1;
       sk[3] = N-2;
       h[0] = 1./12.;
       h[1] =-2./3.;
       h[2] = 2./3.;
       h[3] = -1./12.;
     }
   else
     {
       sk[0] = k-2;
       sk[1] = k-1;
       sk[2] = k+1;
       sk[3] = k+2;
       h[0] = 1./12.;
       h[1] =-2./3.;
       h[2] = 2./3.;
       h[3] =-1./12.;
     }  

   for (i=0; i <P; i++)
     {
      /* if (sk[i] < 0 || sk[i] > N)
         {
           fprintf(stderr,"ERROR: deriv(): sk[%d]=%d\n", i, sk[i]);
           printf("ERRO no calculo da derivada Simétrica\n");
           exit(-1);
         } */
       if (flag==1)
	 v[i] = u[sk[i]][jy]; 
       else
	 v[i] = u[ix][sk[i]];     
     }
   /* compute the derivative */
   du = 0;
   for (i=0; i < P; i++)
     du += v[i] * h[i];

   du /= (1*dxy);

   return du;
 }