#include <stdio.h>
#include <iostream.h>
 
extern "C" {
#include <math.h>
#include "matrix.h"
#include "matrix2.h"
#include "matlab.h"
}

// 	function prototypes, from s3mi.map:
#include "s3mo.inc"


//	global variables:
// The following are in  #include "s3glob.inc" with external linkage
 double a, b, h;
 int m, M, N, j, k;
#undef PI
 double PI = 3.14159265358979323846264338328;
 VEC   *u_vec;
 

// preliminary subroutines:

double sinc(double x)
{
  if (fabs(x) < 0.04) 
    return 1.0-PI*PI*x*x/6+pow(PI,4.0)*pow(x,4.0)/120;
  else
    return sin(PI*x)/(PI*x);
} // ok [Nov-12-1994]


double sinc_p(double x) 	// first derivative
{
  if (fabs(x) < 0.04)
    return -PI*PI*x/3+pow(PI,4.0)*x*x*x/30-pow(PI,6.0)*pow(x,5.0)/840;
  else
    return cos(PI*x)/x-sin(PI*x)/PI/(x*x);
} // ok [Nov-12-1994]


double sinc_pp(double x)	// second derivative
{
  if (fabs(x) < 0.04)
    return -PI*PI/3+pow(PI,4.0)*x*x/10-pow(PI,6.0)*pow(x,4.0)/168;
  else
    return  -sin(PI*x)*PI/x-2.0*cos(PI*x)/(x*x)+2.0*sin(PI*x)/PI/(x*x*x);
} // ok [Nov-12-1994]


//	return a vector pair (x, ux) with plot points from ffrom to tto:
void plot_u(VEC *x, VEC *ux, double ffrom, double tto, int plotpoints)
{
  int i = 0;
  double xx;
  
  for (i=0; i<plotpoints; i++) {
    xx = ffrom + (tto - ffrom)/(plotpoints-1)*i;
    x->ve[i] = xx;
    ux->ve[i] = u(xx);
  }  
}

void plot_u_prime(VEC *x, VEC *upx, double ffrom, double tto, int plotpoints)
{
  int i = 0;
  double xx;
  
  for (i=0; i<plotpoints; i++) {
    xx = ffrom + (tto - ffrom)/(plotpoints-1)*i;
    x->ve[i] = xx;
    upx->ve[i] = u_p(xx);
  }  
}


int main()
{
  int i;
  double x;

//   // 	Simple test:
//   cout << "Test values for sinc(x) and first two derivatives " << endl;
//   for (x = 0.0; x < 6.28; x += 0.1) {
//     cout << x << "  " << sinc(x) << "   " << sinc_p(x) << "   " <<
//       sinc_pp(x) << endl;
//   }
  

  // 	system sizes:
  cout << "Enter M and N: " ;
  cin >> M;
  cin >> N;

  cout << "Enter a and b " ;
  cin >> a;
  cin >> b;

  h = (PI/sqrt(2*M));
  m = M+N+1;

  VEC   *rrhs;
  MAT	*llhs;
  PERM	*pivot;

  // 	init needed arrays
  llhs = m_get(m+2,m+2);	
  pivot = px_get(m+2);
  rrhs = v_get(m+2);

  // 	set up system:
  int indj, indk;
  double xj;
  for (j=-M-1; j <=N+1; j++) {
    indj = j+M+2;
    xj = (psi(j*h));

    //	the 1st spline:
    k = -M-1;
    indk = k+M+2;
    llhs->me[indj-1][indk-1] = coe_spline1(xj); // zero-relative indexing

    // 	the sinc series:
    for (k=-M; k<=N; k++) {
      indk = k+M+2;
      llhs->me[indj-1][indk-1] = coe_sincs(xj); // zero-relative indexing
    }
    // 	the second spline:
    k = N+1;
    indk = k+M+2;
    llhs->me[indj-1][indk-1] = coe_spline2(xj); // zero-relative indexing

    //	The right side:
    rrhs->ve[indj-1] = coe_rhs(xj);
  }
  
  // 	Solve system:
  LUfactor(llhs,pivot);
  u_vec = LUsolve(llhs,pivot,rrhs,VNULL);

  cout << "The coefficient vector:" << endl;
  v_output(u_vec);

  // 	Provide graphs of u and u' :
  {
    VEC *x, *ux;
    double ffrom, tto;
    int plotpoints, i;

    cout << "number of plotpoints?" << endl;
    cin >> plotpoints;
    // 	allocate room:
    x = v_get(plotpoints);
    ux = v_get(plotpoints);

    cout << "Left limit, right limit ? " << endl;
    cin >> ffrom;
    cin >> tto;


    FILE *fp;
    if ( (fp=fopen("tmp.u","w")) == (FILE *) NULL ) {
      printf("Cannot open save file \"%s\"\n", "tmp.u");
      exit(1);
    }
    cout << "x, u(x) for plotting u(x): " << endl;
    plot_u(x,ux, ffrom, tto, plotpoints);
    for (i=0; i<plotpoints; i++)
      fprintf(fp, "  %12.5e  %12.5e \n", x->ve[i], ux->ve[i]);
    fclose(fp);
    

    if ( (fp=fopen("tmp.up","w")) == (FILE *) NULL ) {
      printf("Cannot open save file \"%s\"\n", "tmp.up");
      exit(1);
    }
    cout << endl << "x, u'(x) for plotting u'(x): " << endl;
    plot_u_prime(x,ux, ffrom, tto, plotpoints);
    for (i=0; i<plotpoints; i++)
      fprintf(fp, "  %12.5e  %12.5e \n", x->ve[i], ux->ve[i]);

    fclose(fp);
  }
  return 0;
}