*
Source:    cd_imp_mq.cpp
                  Conveccion - Difusion en 2-d 
                  Metodo implicito 
                  Nucleo Multicuadrico
Compile:
                   c++ -o cdimq cd_imp_mq.cpp -lm 

Example:
                    ./cdimq 30 1 0.1 10 1  0.2 0.001 0.1

Output:
-----------------------------------------------------
Convection-Difussion - Implicit - MQ
n           = 30 , total(nxn) =900
a           = 1.000000
b           = 0.100000
convectivo  = 10.000000
difusivo    = 1.000000
c           = 0.200000
stept       = 0.001000
tmax        = 0.100000

t=0.100000, e_rms       = 0.000314
t=0.100000, e_inf       = 0.001146
-----------------------------------------------------

Los siguientes archivos de datos que
se pueden ver con GNUPLOT
          
            splot  'real'
            splot  'numerica'
            splot  'realnumerica'
            splot  'realt'
            splot  'numericat'
            splot  'realnumericat'

real   -> corresponde a la solucion analitica de la EDP en t=0
realt  -> corresponde a la solucion analitica de la EDP en t=tmax
numerica   -> corresponde a la solucion numerica de la EDP en t=0
numericat   -> corresponde a la solucion numerica de la EDP en t=tmax
realnumerica ->     analitica - numerica  en t=0
realnumericat ->     analitica - numerica  en t=tmax

Autor: Jose Antonio Muñoz Gómez
Supervisión: Pedro González Casanova


*/

#include
#include
#include

#define INTERIOR       0
#define FRONTERA_REAL  1

#include "../libs/jag_mem.cpp"
#include "../libs/jag_phi.cpp"
#include "../libs/jag_vector.cpp"
#include "../libs/jag_cond.cpp"
#include "../libs/qSort.cpp"
#include "../libs/jag_file.cpp"


#include "cd_edp_ut.cpp"

#include "../libs/varios1.cpp"


//--------------------------------------------------------------------
void reconstruye(double *x,double *y, double *lambda,double c,int N, double *sol_numerica)
{
int     i,j;
double xc,yc;
double s;

 for(i=0;i       \n\n");
   exit(1);     
 }
 
    n = atoi(argv[1]);
a_cd = atof(argv[2]);  
b_cd = atof(argv[3]);  
Convectivo = atof(argv[4]);  
Difusivo       = atof(argv[5]);  
c  = atof(argv[6]);  
stept = atof(argv[7]);  
tmax  = atof(argv[8]);  

 tini = 0.0;
 printf("-----------------------------------------------------\n");
 printf("n           = %d , total(nxn) =%d\n",n,n*n);
 printf("a           = %f\n",a_cd);   
 printf("b           = %f\n",b_cd);   
 printf("convectivo  = %f\n",Convectivo);   
 printf("difusivo    = %f\n",Difusivo);   
 printf("c           = %f\n",c);      
 printf("stept       = %f\n",stept);   
 printf("tmax        = %f\n",tmax);    



//Creo el grid cartesiano
 xmin =  0;
 ymin =  0;
 xmax =  1;
 ymax =  1;

  malla  = make_mesh(xmin,xmax,ymin,ymax,n);
  
  N = malla.N;
  x = malla.x;
  y = malla.y;
tipo= malla.tipo;  
 ni = malla.ni; 
 
 
//Reservo la memoria dinamica
  A           = make_dmatrix(N,N);
  M           = make_dmatrix(N,N);
  P           = make_dmatrix(N,N);
  f           = make_dvector(N);
  lambda      = make_dvector(N);
  pivot_index = make_ivector(N);
  
  sol_analitica = make_dvector(malla.N);
  sol_numerica  = make_dvector(malla.N);
  sol_analitica_numerica = make_dvector(malla.N);   
   
  
//Creo la matriz de interpolacion:  condicion inicial
for(int i=0;i=ni){ //frontera
       M[i][j] = u;
     }     
   } 
    
//Creo la matriz-P  
  for(int i=0;i