#include "phead.h"

#define P_max(A,B) ((A)>(B)? (A):(B))

#define LINE printf("*************************************\n")

#ifdef VMS
#include "[-.source]bvalue.c"
#else
#include "../source/bvalue.c"
#endif

/*extern Pfloat bvalue(Pfloat *a, Pfloat* b, Pint C, Pint d, Pint f, Pfloat g)
{
return (Pfloat) 0.0;
}
*/

void outmat(Pchar *cc, Pfloat* bcoef, Pint nx, Pint ny)
{
  Pint ix,iy;
  printf("Dump of matrix %s\n nx:%3ld, ny:%3ld\n",cc,nx,ny);
  for (ix=1; ix<=nx ;ix++){
      for (iy=1; iy<=ny ;iy++) printf("%9.3f",bcoef[ix +(iy-1)*nx -1]);
      printf("\n");
  }
}
  

Pfloat gg(Pfloat x,Pfloat y)
{ Pfloat r1,r2,r3,r4;
    r2 = x - 3.5f;
    r1 = P_max(r2,0.0f);
    r4 = y  - 3.0f;
    r3 = P_max(r4,0.0f);
    return (Pfloat)(r1 * r1) + (r3 * r3 * r3);
}
main(void)
{
    /* Local variables */
    Pfloat taux[7], tauy[6], work1[42], work2[42], work3[42];
    Pint i, j;
    Pfloat bcoef[42]	/* was [7][6] */;
    Pint iflag, mflag, lefty;
    Pint jj, kx, ky, nx, ny;
    Pfloat tx[10], ty[10];
    Pfloat xx,yy;
    
    nx = 7;
    kx = 3;
    ny = 6;
    ky = 4;
    
    for (i = 1; i <= nx; i++) {
	taux[i - 1] = (Pfloat) i;
    }
    for (i = 1; i <= kx; i++) {
	tx[i - 1] = taux[0];
	tx[nx + i - 1] = taux[nx - 1];
    }

    for (i = kx + 1; i <= nx; i++) {
	tx[i - 1] = (taux[i - kx] + taux[i - kx + 1]) / 2.0f;
    }

    for (i = 1; i <= ny; i++) {
	tauy[i - 1] = (Pfloat) i;
    }

    for (i = 1; i <= ky; i++) {
	ty[i - 1] = tauy[0];
	ty[ny + i - 1] = tauy[ny - 1];
    }

    for (i = ky + 1; i <= ny; i++) {
	ty[i - 1] = tauy[i - ky + 1];
    }

    printf("Given Data\n\n          ");
    for (i = 1; i <= ny; i++) printf("%10.1f",tauy[i-1]);
    printf("\n");

    for (i = 1; i <= nx; i++) {
	printf("%10.1f",taux[i-1]);
	for (j = 1; j <= ny; j++) {
	    bcoef[i + j * 7 - 8] = gg(taux[i-1],tauy[j-1]);
	    printf("%10.3f",bcoef[i + j * 7 - 8]);
	}
	printf("\n");
    }

    fflush(stdout);

    spli2d(taux, bcoef, tx, nx, kx, 1, 1, ny, work2, work3, work1, &iflag);
    outmat("Work1",work1,nx,ny);
    spli2d(tauy, work1, ty, ny, ky, 1, 1, nx, work2, work3, bcoef, &iflag);
    outmat("\n\nbcoef",bcoef,nx,ny);
    
    printf("\n");
    printf("\n");

    for (j = 1; j <= ny; j++) {
	interv(ty, ny, tauy[j - 1], &lefty, &mflag);
	for (i = 1; i <= nx; i++) {
	    for (jj = 1; jj <= ky; jj++) {
		work2[jj - 1] = bvalue(tx, &bcoef[(lefty - ky + jj) * 7 - 7],
			 nx, kx, taux[i - 1]);
	    }
	    yy = gg(taux[i-1],tauy[j-1]);
	    xx = bvalue(&ty[lefty - ky], work2, ny, ky, tauy[j - 1]);
	    work1[i + j * 7 - 8] = yy - xx;
	}
    }

    printf("Interpolation Error\n\n          ");
    for (i = 1; i <= ny; i++) printf("%10.3f",tauy[i-1]);
    printf("\n");

    for (i = 1; i <= nx; i++) {
        printf("%10.1f",taux[i-1]);
	for (j = 1; j <= ny; j++) printf("%10.3f",work1[i + j * 7 - 8] );
    	printf("\n");
    }
     
}

