//Hi folks! i need to know hot to declare an array
//
//I know how to dogie


/*
    ##############################################################################################################
        FUNCIÓN PARA CALCULAR LA DERIVADA PRIMERA DE UNA FUNCIÓN USANDO LAS FÓRMULAS VISTAS EN CLASE
          - LOS PUNTOS VIENEN DADOS EN EL VECTOR x.
          - LA FUNCIÓN f(x) VIENE DADA POR EL VECTOR fx,
          - HAY QUE TENER EN CUENTA LA GESTIÓN DE LOS BORDES, ES DECIR, DE LOS VALORES DE LOS EXTREMOS DEL VECTOR:
          - BORDES: CALCULAR DERIVADA CON 2 PUNTOS. ES DECIR : f'(xi)=(fx(xj)-fx(xi))/(xj-xi)
          - RESTO DE ELEMENTOS: CALCULAR DERIVADA CON 3 PUNTOS. ES DECIR:
          f'(xi)=((xi-xl)(f(xr)-f(xi))/(xr-xi) + (xr-xi)(f(xi)-f(xl))/(xi-xl))  / (xr-xl)
          - LA DERIVADA SE DEVUELVE EN EL VECTOR DE ENTRADA/SALIDA fd.
          - LA FUNCIÓN DEVUELVE 0 SI TERMINA CORRECTAMENTE, Y -1 EN CASO DE ERROR.
    ################################################################################################################
  */

int mn_derivada(Array1D <real> x, Array1D <real> fx, Array1D <real> &fd) {
  int Npuntos = x.dim();
  //HACER ALUMNO
  int Fpuntos = fx.dim();
  if(Npuntos != Fpuntos) {
    return -1;
  } else {
    for (int i = 0; i < Npuntos; i++) {
      if(i == 0) {
        fd[i] = (fx[i+1] - fx[i])/(x[i+1] - x[i]);
      } else if(i == Npuntos-1) {
        fd[i] = (fx[i-1] - fx[i])/(x[i-1] - x[i]);
      } else {
        fd[i]=((((x[i] - x[i-1]) * (fx[i+1] - fx[i])) / (x[i+1] - x[i])) + (((x[i+1] + x[i]) * (fx[i] - fx[i-1])) / (x[i] - x[i-1]))) / (x[i+1] - x[i-1]);
      }
    }
  }
  return 0;
}



// #############################################################
// FUNCIÓN PARA CALCULAR UNA INTEGRAL POR EL MÉTODO DEL TRAPECIO
// UTILIZANDO N INTERVALOS EN UN INTERVALO [a,b]
// #############################################################

real mn_trapecio(
  real (*f)(const real x),  // Función que se integra
  const real a,       // Extremo izquierdo del intervalo
  const real b,       // Extremo derecho del intervalo
  const int N)        // Número de intervalos para calcular la integral
{
  //HACER ALUMNO
  real integral = 0;                                  // Variable donde se acumula el valor de la integral
  real funcion = (*f)(a);                               // Evaluación de la función en a
  real h = (b-a)/(2.*N);                                // Factor multiplicativo
  for (int i = 0; i < N; i++) {                           // Bucle para acumular el valor de la integral
    real ext_dch = a + (i+1) * (b-a)/N;                       // Extremo derecho subintervalo
    real ev_funcion = (*f)(ext_dch);                        // Evaluación de la función en el extremo derecho
    integral += h * (funcion + ev_funcion);                     // Acumulación valor integral
    funcion = ev_funcion;                             // Actualización valor integral extremo izquierdo
  }
  return integral;      // Se devuelve el valor de la integral
}

// ###########################################################################
// FUNCION PARA CALCULAR UNA INTEGRAL POR ITERACIONES DEL METODO DEL TRAPECIO
// SE PARTE DE N=1 INTERVALOS Y SE CALCULA LA INTEGRAL CON LA FUNCION ANTERIOR
// SE ACTUALIZA N MULTIPLICANDOLO POR 2 Y SE CALCULA DE NUEVO LA INTEGRAL
// SE REALIZA ESTE PROCESO HASTA QUE LA DISTANCIA RELATIVA ENTRE LOS 2
// VALORES DE LA INTEGRAL ES MENOR QUE LA TOLERANCIA. EN Niter  SE DEVUELVE
// EL NUMERO DE ITERACIONES QUE SE REALIZA DE ESTE PROCESO
// ###########################################################################
mn_trapecio(
  real (*f)(const real x),  // Función que se integra
  const real a,       // Extremo izquierda del intervalo
  const real b,       // Extremo derecho del intervalo
  int &Niter,         // Variable de salida con el número de iteraciones realizadas
  const real TOL)       // Tolerancia para controlar la convergencia
{
  // HACER ALUMNO
  int N = 1;
  real integral = mn_trapecio((*f),a,b,N);
  real error = TOL + 1;
  Niter = 0;
  while (error > TOL) {
    if(N > 1e7) break;    // Fijamos "a priori" el número máximo de intervalos a 1e7
    Niter++;
    N *= 2;
    real integral2 = mn_trapecio((*f),a,b,N);
    error = mn_distancia(integral, integral2);
    integral = integral2;
  }
  return (integral);
}



// ###########################################################
// FUNCIÓN PARA CALCULAR UNA INTEGRAL POR EL MÉTODO DE SIMPSON
// ###########################################################
real mn_simpson(
  real (*f)(const real x),  // Función que se integral
  const real a,       // Extremo izquierdo del intervalo
  const real b,       // Extremo derecho del intervalo
  const int N)        // Número de intervalos para calcular la integral
{
  //HACER ALUMNO
  real integral = 0;                                  // Variable donde se acumula el valor de la integral
  real funcion = (*f)(a);                               // Evaluación de la función en a
  real h = (b-a)/(6.*N);                                // Factor multiplicativo
  for (int i = 0; i < N; i++) {                           // Bucle para acumular el valor de la integral
    real ext_izq = a + i * (b-a)/N;                         // Extremo izquierdo subintervalo
    real ext_dch = a + (i+1) * (b-a)/N;                       // Extremo derecho subintervalo
    real ev_funcion = (*f)(ext_dch);                        // Evaluación de la función en el extremo derecho
    integral += h * (funcion + ev_funcion + 4 * (*f)((ext_izq + ext_dch)/2.));    // Acumulación valor integral
    funcion = ev_funcion;                             // Actualización valor integral extremo izquierdo
  }
  return integral;      // Se devuelve el valor de la integral
}

// ########################################################################
// FUNCIÓN PARA CALCULAR UNA INTEGRAL POR ITERACIONES DEL MÉTODO DE SIMPSON
// ########################################################################
real mn_simpson(
  real (*f)(const real x),  // Función que se integra
  const real a,       // Extremo izquierda del intervalo
  const real b,       // Extremo derecho del intervalo
  int &Niter,         // Variable de salida con el número de iteraciones realizadas
  const real TOL)       // Tolerancia para controlar la convergencia
{
  // HACER ALUMNO
  int N = 1;
  real integral = mn_simpson((*f),a,b,N);
  real error = TOL + 1;
  Niter = 0;
  while (error > TOL) {
    if(N > 1e7) break;    // Fijamos "a priori" el número máximo de intervalos a 1e7
    Niter++;
    N *= 2;
    real integral2 = mn_simpson((*f),a,b,N);
    error = mn_distancia(integral, integral2);
    integral = integral2;
  }
  return (integral);
}



// ###############################################################
// FUNCIÓN PARA CALCULAR UNA INTEGRAL POR EL MÉTODO DEL RECTÁNGULO
// ###############################################################
real mn_rectangulo(
  real (*f)(const real x),  // Función que se integral
  const real a,       // Extremo izquierdo del intervalo
  const real b,       // Extremo derecho del intervalo
  const int N)        // Número de intervalos para calcular la integral
{
  //HACER ALUMNO
  real integral = 0;                                  // Variable donde se acumula el valor de la integral
  real funcion = (*f)(a);                               // Evaluación de la función en a
  real h = (b-a);                                   // Factor multiplicativo
  for (int i = 0; i < N; i++) {                           // Bucle para acumular el valor de la integral
    real ext_izq = a + i * (b-a)/N;                         // Extremo izquierdo subintervalo
    real ext_dch = a + (i+1) * (b-a)/N;                       // Extremo derecho subintervalo
    real ev_funcion = (*f)(ext_dch);                        // Evaluación de la función en el extremo derecho
    integral += h * ((*f)((ext_izq + ext_dch)/2.));                 // Acumulación valor integral
    funcion = ev_funcion;                             // Actualización valor integral extremo izquierdo
  }
  return integral;      // Se devuelve el valor de la integral
}



// ############################################################################
// FUNCIÓN PARA CALCULAR UNA INTEGRAL POR ITERACIONES DEL MÉTODO DEL RECTÁNGULO
// ############################################################################
real mn_rectangulo(
  real (*f)(const real x),  // Función que se integra
  const real a,       // Extremo izquierda del intervalo
  const real b,       // Extremo derecho del intervalo
  int &Niter,         // Variable de salida con el número de iteraciones realizadas
  const real TOL)       // Tolerancia para controlar la convergencia
{
  // HACER ALUMNO
  int N = 1;
  real integral = mn_rectangulo((*f),a,b,N);
  real error = TOL + 1;
  Niter = 0;
  while (error > TOL) {
    if(N > 1e7) break;    // Fijamos "a priori" el número máximo de intervalos a 1e7
    Niter++;
    N *= 2;
    real integral2 = mn_rectangulo((*f),a,b,N);
    error = mn_distancia(integral, integral2);
    integral = integral2;
  }
  return (integral);
}

/************************************************************
FUNCION PARA CALCULAR UNA INTEGRAL POR EL METODO DE CUADRATURA DE GAUSS UTILIZANDO
UN VECTOR DE PUNTOS Y PESOS PARA EL INTERVALO [-1,1]. LOS VALORES DE PUNTOS Y PESOS
VIENEN EN LAS MATRICES w Y x, POR EJEMPLO PARA LA FÓRMULA DE 2 PUNTOS LOS PESOS
SERÍAN w[1][0] Y w[1][1] y LOS PUNTOS SERÍAN x[1][0] Y x[1][1]

NOTA : HAY QUE TENER EN CUENTA EL CAMBIO DE VARIABLE VISTO EN CLASE PARA PASAR UNA
INTEGRAL EN UN INTERVALO [a,b] AL INTERVALO [-1,1]

EN ESTE CASO LA INTEGRAL ES EN 2 DIMENSIONES Y POR TANTO LA FUNCIÓN A INTEGRAR DEPENDE
DE 2 VARIABLES Y HAY QUE AJUSTAR LA INTEGRAL TAL Y COMO SE EXPLICÓ EN CLASE. EL Nº DE
PUNTOS DE CUADRATURA A UTILIZAR ES EL MISMO EN LA VARIABLE x E y Y VIENE DADO EL PARÁMETRO
DE ENTRADA N
************************************************************/

real mn_cuadratura(
real (*f)(real x,real y) /* funcion que se integra */,
real a /* extremo izquierdo del intervalo para x*/,
real b /* extremo derecho del invervalo para x*/,
real c /* extremo inferior del intervalo para y*/,
real d /* extremo superior del invervalo para y*/,
Array2D< real > &x, /* matriz con los puntos de la tabla de cuadratura para el intervalo [-1,1]*/
Array2D< real > &w, /* matriz con los pesos de la tabla de cuadratura para el intervalo [-1,1] */
int N) /* Nº de puntos a usar en la fórmula de cuadratura */
{
   // HACER ALUMNO

   real suma = 0.; // variable donde guardare el resultado
   real paso = 0.; // variables auxiliares que creo para no tener una linea demasiado larga
   real paso2 = 0.; // a la hora de hacer los calculos
   real imagen = 0.;
    for(int k = 0; k<N; k++){ // bucle de la primera integral
        for(int j = 0; j<N; j++){ //bucle de la segunda integral
            paso = ((b-a)*x[N-1][k]+b+a)/2.; // valor que depende del primer intervalo [a,b]
            paso2 = ((d-c)*x[N-1][j]+d+c)/2.; // valor que depende del segundo intervalo [c,d]
            imagen = (*f)(paso,paso2); // calculo la funcion en esos puntos
            suma += w[N-1][k]*w[N-1][j]*((b-a)/2.)*((d-c)/2.)*imagen; // calculo el sumatorio
        }
   }
   return suma;
}

/************************************************************
FUNCION PARA CALCULAR UNA INTEGRAL POR EL METODO DE CUADRATURA DE GAUSS UTILIZANDO
UN VECTOR DE PUNTOS Y PESOS PARA EL INTERVALO [-1,1]

NOTA : HAY QUE TENER EN CUENTA EL CAMBIO DE VARIABLE VISTO EN CLASE PARA PASAR UNA
INTEGRAL EN UN INTERVALO [a,b] AL INTERVALO [-1,1]

LA FUNCIÓN DEVUELVE EL VALOR DE LA INTEGRAL.
************************************************************/
real mn_cuadratura(
real (*f)(real x) /* funcion que se integra */,
real a /* extremo izquierdo del intervalo */,
real b /* extremo derecho del invervalo*/,
Array1D< real > &x, /* vector con los puntos de la tabla de cuadratura para el intervalo [-1,1]*/
Array1D< real > &w) /* vector con los pesos de la tabla de cuadratura para el intervalo [-1,1] */

{
  //A IMPLEMENTAR POR EL ALUMNO.
  real suma = 0.;
  real imagen = 0.;
  real paso = 0.;
  int N = x.dim();

  for (int k = 0; k<N; k++) {
      paso = ((b-a)*x[k]+b+a)/2.;
      imagen = (*f)(paso);
      suma = suma + w[k]*((b-a)/2.)*imagen;
  }

  return suma;
}


//Funcion que realiza la aproximacion de una derivada con 3 ptos
real mn_derivada3(real (*f)(real x),real xr, real xi){

real h = xr-xi;
real derivada;

derivada =(f(xi+h)-f(xi-h))/2;
return derivada;

}

// Funcion que realiza la derivada con 2 ptos.
real mn_deriva2(real (*f)(real x),real xr, real xi){

if (xr == xi) {return -1.;}
real derivada = (f(xr)-f(xi))/(xr-xi);
return derivada;
}

//============================================================================================
// FUNCION QUE CALCULA UNA INTEGRAL POR ITERACIONES DEL METODO DE SIMPSON:
//  - SE PARTE DE N=1 INTERVALOS Y SE CALCULA LA INTEGRAL CON LA FUNCION ANTERIOR (mn_simpson)
//  - SE ACTUALIZA N MULTIPLICANDOLO POR 2 Y SE CALCULA DE NUEVO LA INTEGRAL
//  - SE REALIZA ESTE PROCESO HASTA QUE LA DISTANCIA RELATIVA ENTRE LOS 2 VALORES DE LA INTEGRAL
//    SEA MENOR QUE LA TOLERANCIA O SE ALCANCE EL NÚMERO MÁXIMO DE INTERVALOS.
//  - EN Niter SE ALMACENA EL NUMERO DE ITERACIONES QUE SE REALIZAN EN ESTE PROCESO.
//============================================================================================
real mn_simpson(
real (*f)(real x) /* funcion que se integra */,
real a /* extremo izquierdo del intervalo */,
real b /* extremo derecho del invervalo*/,
real TOL /* tolerancia para controlar la convergencia */,
real Imax /* número máximo de intervalos */,
int &Niter /* variable de salida con el numero de iteraciones realizadas */
)
{
  /* A IMPLEMENTAR POR EL ALUMNO */
  //NOTA: CALCULAR EL ERROR CON mn_distancia()
  real integral = 0.;
  real integral1 = 0.;
  Niter=0;
  real x= (b-a)/Imax;
  for (int j=1;j<Imax;j++){
     integral= mn_simpson(f,a,b,j);
     integral1= mn_simpson(f,a,b,(2*j));
     Niter++;
     if (mn_distancia(integral,integral1)<TOL || j==Imax-1){
        break;
     }
  }
}