/*
 *  Riccflowcode.c
 *  
 *
 *  Created by Christopher Herzog on 11/24/06.
 *  Copyright 2006. All rights reserved.
 *
 */


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include "Ricciflowcode.h"

int main(int argc, char *argv)
{
extern char *optarg;
extern double DELTA;
extern int NMAX;
/*extern double curvol;*/
extern double FLOW_PARAM;

double **tf;
double **newfunct;

FILE *fpin, *fpout;

fpin = NULL;
fpout = NULL;

int i,c,j, myflag, counter, pcounter;
int myprint, Ntest;
double WIDTH, MAX_LOOPS, fshift, error, DIFF_ERROR;

double foo, discerror;
int bar;

double oldshift;


WIDTH = WIDTH_DEFAULT;
NMAX = NMAX_DEFAULT;
MAX_LOOPS = MAX_LOOPS_DEFAULT;
myprint = MY_PRINT_DEFAULT;

/*curvol = 0.5;*/
FLOW_PARAM = FLOW_PARAM_DEFAULT;


/* l controls number of loops
 * w controls grid width
 * N controls grid size 
 * i indicates input grid
 * o indicates output file for grid
 * d indicates how many iterations to run before displaying run data
 * R controls Ricci flow time steps
 */

	while ((c = getopt(argc, argv, "l:w:N:i:o:d:R:")) != -1){
		switch (c) {
		case 'l':
			MAX_LOOPS = atoi(optarg);
			break;
		case 'w':
			WIDTH = atof(optarg);
			break;
		case 'N':
			NMAX = atoi(optarg);
			break;
		case 'i':
			fpin = fopen(optarg, "r");
			if(fpin==NULL){
				printf("Couldn't open input file: exiting\n");
				exit(1);
			}
			break;
		case 'o':
			fpout = fopen(optarg, "w");
			if(fpout==NULL){
				printf("Couldn't open output file: exiting\n");
				exit(1);
			}
			break;
		case 'd':
			myprint = atoi(optarg);
			/*printf("%d \n", myprint);*/
			break;
		/*case 'v':
			curvol = atof(optarg);
			break;*/
		case 'R':
			FLOW_PARAM = atof(optarg);
			break;
		default:
			printf("Illegal Option\n");
			exit(1);
			break;
		}		
	}

DELTA = WIDTH / (NMAX-1);
DIFF_ERROR = DIFF_ERROR_RAW * pow(NMAX,2);

/*printf("%f \n", Fpottest(0.7,0.6));
 *exit(1);
 */

if(fpout == NULL) fpout = fopen(DEFAULT_OUT, "w");

if(fpout == NULL) {
	printf("Couldn't open default output file: exiting\n");
	exit(1);
}

if(fpout == fpin){
	printf("Input file is output file: exiting \n");
	exit(1);
}

tf = calloc(NMAX+2, sizeof(double *));
newfunct = calloc(NMAX+2, sizeof(double *));
for(i=0; i<= NMAX+1; i++){
 /* printf("%d \n", i);*/
 tf[i] = (double *) calloc(NMAX+2, sizeof(double));
 newfunct[i] = (double *) calloc(NMAX+2, sizeof(double));
 if(tf[i]==NULL){
	printf("Memory allocation for grid failed. \n");
	exit(1);
 }
}

if(fpin == NULL) {
	printf("Using canonical potential as initial value. \n");
	for(i=1; i<=NMAX+1; i++){
		for(j=1; j<=NMAX+1; j++){
			/*printf("%d, %d \n", i,j);*/
			tf[i][j] = Fpottest(xicoord(i), etacoord(j));
		}
	}
	boundaryupdate(tf);
}
else {
	printf("Using input file. \n");
	fscanf(fpin, "%d\n",&bar); 
	fscanf(fpin, "%f\n",&foo);
	
	fscanf(fpin, "array = {");
	/*if(Ntest != NMAX){
		printf("input grid of size %d does not match nmax=%d: exiting \n", Ntest, NMAX);
		exit(1);
	}*/
	for(i=0; i<=NMAX+1; i++){
		fscanf(fpin, "{");
		c=0;
		for(j=0; j<=NMAX+1; j++){
			if(j==NMAX+1) fscanf(fpin, "%lf", &tf[i][j]);
			else fscanf(fpin, "%lf, ", &tf[i][j]);
			c++;
			if(c==5){
				fscanf(fpin, "\n");
				c=0;
			}
		}
		if(i==NMAX+1) fscanf(fpin, "}");
		else fscanf(fpin, "},\n");
	}
	fscanf(fpin, "}");
	fclose(fpin);
}

printf("NMAX %d : DELTA %f \n", NMAX, DELTA);

myflag = 1;
counter = 1;
pcounter = 1;

/*printf("%d, %d \n", pcounter, printoutvar);*/

oldshift = 0;
fshift = 0;

while(myflag){
	boundaryupdate(tf);
	oldshift = fshift;
	fshift = Fijupdate(tf, newfunct, &oldshift);
	error = 0;
	for(i=1; i<=NMAX; i++){
		for(j=1; j<=NMAX; j++){
			error += fabs(tf[i][j]-newfunct[i][j]);
		}
	}
	if(error < DIFF_ERROR || counter > MAX_LOOPS){
		myflag = 0;
	}
	if(pcounter == myprint || counter == 2){
		discerror = calcdiscerror(tf,fshift);
		printf("%d : %e : %e : %e : %-16.14e \n", counter, error, fshift, discerror, 
					0.5-volcomp(tf, fshift));
		if(pcounter == myprint) pcounter = 0;
	}
	/*following loop for Jacobi differencing only */
	/*for(i=0; i<=NMAX+1; i++){
		for(j=0; j<=NMAX+1; j++){
			tf[i][j] = newfunct[i][j];
		}
	}*/
	counter++;
	pcounter++;
}

boundaryupdate(tf);

fprintf(fpout,"%d\n",NMAX);
fprintf(fpout,"%f\n",DELTA);

fprintf(fpout, "array = {");
for(i=0; i<=NMAX+1; i++){
	fprintf(fpout, "{");
	c=0;
	for(j=0; j<=NMAX+1; j++){
		if(j==NMAX+1) fprintf(fpout, "%-16.14f", tf[i][j]-fshift/2.);
		else fprintf(fpout, "%-16.14f, ", tf[i][j]-fshift/2.);
		c++;
		if(c==5){
			fprintf(fpout, "\n");
			c=0;
		}
	}
	if(i==NMAX+1) fprintf(fpout, "}");
	else fprintf(fpout, "},\n");
}
fprintf(fpout, "}");
fclose(fpout);

printf("%d : %e : %-16.14f: %-16.14f \n", counter-1, error, fshift, volcomp(tf,fshift));

}

double calcdiscerror(double **tf, double oldshift){
	extern int NMAX;
	int i,j,flag;
	double error;
	
	error=0;
	for(i=2; i<=NMAX; i++){
		for(j=2; j<=NMAX; j++){
		
		flag = 0;

		error += fabs(discfunct(i,j, flag, tf)-oldshift);
		
		
		}
	}

	return(error/ (double) pow(NMAX,2));

}

double volcomp(double **tf, double fshift){

	int i,j;
	double cadj, newvar;
	
	cadj = 0;

	for(i=1; i<=NMAX; i++){
		for(j=1; j<=NMAX; j++){

		newvar = tf[i][j];
		if(xicoord(i)>0 && xicoord(i) <= 1.0- DELTA/2.0 &&
			  etacoord(j)>0 && etacoord(j) <= 1.0 - DELTA/2.0){
						cadj += exp(-2*newvar + fshift)*xicoord(i)*etacoord(j);
			}
			if(xicoord(i) <=1.0+DELTA/2.0 && xicoord(i) > 1.0-DELTA/2.0 && etacoord(j) <= 1.0-DELTA/2.){
				cadj += (-xicoord(i)+DELTA/2.0+1.0)/DELTA*exp(-2*newvar+fshift)*xicoord(i)*etacoord(j);
			}
			if(etacoord(j) <=1.0+DELTA/2 && etacoord(j) > 1.0-DELTA/2.0 && xicoord(i) <= 1.0-DELTA/2.){
				cadj += (-etacoord(j)+DELTA/2+1.0)/DELTA*exp(-2*newvar+fshift)*xicoord(i)*etacoord(j);
			}
			if(xicoord(i)<=1.0+DELTA/2.0 && xicoord(i) > 1.0-DELTA/2.0 &&
					etacoord(j) <=1.0+DELTA/2.0 && etacoord(j) > 1.0-DELTA/2.0 ){
				cadj += pow((-etacoord(j)+DELTA/2+1.0)/DELTA,2)*exp(-2*newvar+fshift)*xicoord(i)*etacoord(j);
			}
		}
	}
	
	return(cadj * pow(DELTA,2));

}

double xicoord(int i){
	extern double DELTA;
	return((i-1)*DELTA);
}


double etacoord(int j){
	extern double DELTA;
	return((j-1)*DELTA);
}

double Fpottest(double xi, double eta){

  double x1, x2, x1new, x2new, error1, error2, delta1, delta2, denom;
  double xi2, eta2;
  int iterations;
  iterations = 0;
  x1 = 0.5; x2 = -0.5;
  delta1=1; delta2=1;
  error1=1; error2=1;
  xi2 = pow(xi,2);
  eta2 = pow(eta,2);	
	
  while((delta1+delta2 > ROOT_ERROR || error1+error2 > ROOT_ERROR) && iterations<50)
  {
	denom = (2*(-1+x1+x2)*(1+x1+x2)*(-3+2*(pow(x1,2)+x1*x2+pow(x2,2))));
	x1new = x1 +((1+x1)*((1-x1)*(-1+x1+x2)*(-eta2*pow(-1+x2,2)*(-1+x1+x2)+(1+x1+x2)*
			(-3+pow(x1,2)+2*x1*x2+3*pow(x2,2)))+(1+x1)*pow(1+x1+x2,2)*
			(-2+pow(x1,2)+2*x1*x2+2*pow(x2,2))*xi2))/denom;
	x2new = x2 -((-1+x2)*(eta2*(-1+x2)*pow(-1+x1+x2,2)*(-2+2*pow(x1,2)+
			2*x1*x2+pow(x2,2))+
			(1+x2)*(-3+6*pow(x1,2)-3*pow(x1,4)+8*x1*x2-8*pow(x1,3)*x2+4*pow(x2,2)-
			8*pow(x1,2)*pow(x2,2)-4*x1*pow(x2,3)-pow(x2,4)+
			pow(1+x1,2)*pow(1+x1+x2,2)*xi2)))/denom;

   
	error1 = fabs((1-x1new)*(1-x1new-x2new)/(1+x1new)/(1+x1new+x2new) - xi2);
	error2 = fabs((1+x2new)*(1+x1new+x2new)/(1-x2new)/(1-x1new-x2new) - eta2);
	delta1 = fabs(x1-x1new);
	delta2 = fabs(x2-x2new);
	/*printf("%f ; %f ; %f ; %f \n", error1, error2, delta1, delta2);*/
	x1 = x1new; x2 = x2new; 
	iterations++;
  }

	return(-1/2.0*(2.0*log(1+x1new)+log(1+x1new+x2new)+
			log(1-x1new-x2new)+2.0*log(1-x2new)));
}

void boundaryupdate(double **tf){
	extern int NMAX;
	int i,j;
	double temp;
	
	/*printf("reached boundary update \n");*/
	
	for(i=1; i<=NMAX; i++){
		tf[0][i] = tf[2][i];
		tf[i][0] = tf[i][2];
		if(xicoord(i)*etacoord(NMAX+1) > 1.0 && xicoord(i) < 1.0 /*(NMAX-1)*DELTA*/){
		  tf[i][NMAX+1] = planarinterp(1.0/etacoord(NMAX+1)/xicoord(i), xicoord(i), tf)+
						log(xicoord(i)*pow(etacoord(NMAX+1),2));
	      tf[NMAX+1][i] = planarinterp(etacoord(i), 1.0/xicoord(NMAX+1)/etacoord(i), tf)+
						log(etacoord(i)*pow(xicoord(NMAX+1),2));
		} else if(xicoord(i)*etacoord(NMAX+1) <= 1.0) {
		  tf[i][NMAX+1] = planarinterp(1.0/etacoord(NMAX+1), xicoord(i)*etacoord(NMAX+1),tf) +
						log(etacoord(NMAX+1));
		  
		  tf[NMAX+1][i] = planarinterp(etacoord(i)*xicoord(NMAX+1), 1.0/xicoord(NMAX+1), tf) +
						log(xicoord(NMAX+1));
	
		} else {
			tf[i][NMAX+1] =  planarinterp(1.0/xicoord(i), 1.0/etacoord(NMAX+1), tf)+
						log(pow(xicoord(i),2)*pow(etacoord(NMAX+1),2));
			
			tf[NMAX+1][i] = planarinterp(1.0/xicoord(NMAX+1), 1.0/etacoord(i), tf)+
						log(pow(xicoord(NMAX+1),2)*pow(etacoord(i),2));
		
		}
		/*printf("%d ", i);*/
	}
	tf[0][0] = tf[2][2];
	/*
	for(i=1; i<=NMAX; i++){
		tf[NMAX+1][i] = tf[NMAX-1][i];
		tf[i][NMAX+1] = tf[i][NMAX-1];
	}*/
	tf[NMAX+1][NMAX+1] = planarinterp(1.0/xicoord(NMAX+1), 1.0/etacoord(NMAX+1), tf)+
						log(pow(xicoord(NMAX+1),2)*pow(etacoord(NMAX+1),2));
	tf[NMAX+1][0] = tf[NMAX+1][2];
	tf[0][NMAX+1] = tf[2][NMAX+1];

}

double planarinterp(double x, double y, double **tf){
	extern double DELTA;
	extern int NMAX;
	int i,j;
	double xmin, ymin, myresult;
	myresult = 0;
	
	double yy[5], yy1[5], yy2[5], yy12[5];

	double ansy, ansy1, ansy2;
	double DELTA2 = DELTA*DELTA;
	
	/*The ROOR_ERROR is because for a few points the
	 *corner might otherwise be at (NMAX+1,NMAX+1).  Confused
	 *about this now...
	 */
	i = floor(x/DELTA + 1 /*- ROOT_ERROR*/);
	j = floor(y/DELTA + 1 /*- ROOT_ERROR*/);
	if(i==NMAX+1) i=NMAX;
	if(j==NMAX+1) j=NMAX;
	xmin = (i-1)*DELTA;
	ymin = (j-1)*DELTA;
	
	/*printf("%d, %d \n", i,j);*/
	
	if(i+1>NMAX+1 || j+1> NMAX+1){
	   printf("Out of bounds in planarinterp: exiting \n");
	   exit(1);
	  }
	
	if(x<xmin -ROOT_ERROR || y < ymin -ROOT_ERROR || x > xmin + DELTA +ROOT_ERROR || y > ymin + DELTA +ROOT_ERROR){
	   printf("x %f or y %f not in square (%f, %f) \n", x, y, xmin, ymin);
	   
	  }
	
	if(i-1<0 || j-1<0 || i+2 > NMAX+1 || j+2 > NMAX+1){ 
		return(tf[i][j]*(1.0+(-x+xmin)/DELTA)*(1.0+(-y+ymin)/DELTA) + tf[i+1][j]*(x-xmin)/DELTA*(1.0+(-y+ymin)/DELTA)+
			tf[i][j+1]*(1.0+(-x+xmin)/DELTA)*(y-ymin)/DELTA + tf[i+1][j+1]*(x-xmin)*(y-ymin)/DELTA/DELTA);
	} else {
	
	yy[1] = tf[i][j];
	yy[4] = tf[i][j+1];
	yy[2] = tf[i+1][j];
	yy[3] = tf[i+1][j+1];

	yy2[1] = (tf[i][j+1] - tf[i][j-1])/2.0/DELTA;
	yy2[4] = (tf[i][j+2] - tf[i][j])/2.0/DELTA;
	yy2[2] = (tf[i+1][j+1] - tf[i+1][j-1])/2.0/DELTA;
	yy2[3] = (tf[i+1][j+2] - tf[i+1][j])/2.0/DELTA;

	yy1[1] = (tf[i+1][j] - tf[i-1][j])/2.0/DELTA;
	yy1[4] = (tf[i+1][j+1] - tf[i-1][j+1])/2.0/DELTA;
	yy1[2] = (tf[i+2][j] - tf[i][j])/2.0/DELTA;
	yy1[3] = (tf[i+2][j+1] - tf[i][j+1])/2.0/DELTA;

	yy12[1] = 1/4.0*(tf[i+1][j+1] - tf[i+1][j-1] - tf[i-1][j+1] + tf[i-1][j-1])/DELTA2;
	yy12[4] = 1/4.0*(tf[i+1][j+2] - tf[i+1][j] - tf[i-1][j+2] + tf[i-1][j])/DELTA2;
	yy12[2] = 1/4.0*(tf[i+2][j+1] - tf[i][j+1] - tf[i+2][j-1] + tf[i][j-1])/DELTA2;
	yy12[3] = 1/4.0*(tf[i+2][j+2] - tf[i][j+2] - tf[i+2][j] + tf[i][j])/DELTA2;

	bcuint(yy,yy1,yy2,yy12,xmin,xmin+DELTA,ymin,ymin+DELTA, x, y, &ansy, &ansy1, &ansy2);
	
	return(ansy);
	
	}

/*	
	if(y > x+ymin-xmin){
	  myresult = (tf[i+1][j+1]-tf[i][j+1])*x/DELTA +
				(tf[i][j+1]-tf[i][j])*y/DELTA +
				(tf[i+1][j+1]*(1-i)+tf[i][j+1]*(i-j)+tf[i][j]*j);
	  return(myresult);
	} else {
	  myresult = (-tf[i][j]+tf[i+1][j])*x/DELTA +
				(-tf[i+1][j]+tf[i+1][j+1])*y/DELTA +
				(tf[i+1][j+1]*(1-j)+tf[i][j]*i + tf[i+1][j]*(j-i));
	  return(myresult);
	}
*/	
}

/* returns the new lattice in tf, the old lattice in newfunct */

double Fijupdate(double **tf, double **newfunct, double *oldshift){
	extern int NMAX;
	extern double DELTA;
	/*extern double curvol;*/
	extern double FLOW_PARAM;
	
	int i,j, flag, iterations;
	int oneint;
	
	double /*xderiv, yderiv,*/ tempvar, newvar, error, epsilon; /*, cadj;*/
	double fshift, temp;
	double oneoverNMAX2 = FLOW_PARAM/pow(NMAX,2);
	
	/*cadj = 0;*/
	
	/*oneint = (int) 1./DELTA + 1;*/
	/*printf("%d \n", oneint);*/
	
	/*if(oneint == NMAX+1) oneint--;*/
	
	for(i=1; i<=NMAX; i++){
		for(j=1; j<=NMAX; j++){
			epsilon = 1;
			flag = 0; iterations = 0;
			if(i==1 && j==1) flag =1;
			if(j==1 && i!=1) flag =2;
			if(i==1 && j!=1) flag =3;
			tempvar = discfunct(i,j,flag,tf);
			if(i == NMAX/2 && j == NMAX/2){
				fshift = tempvar;
			}
			newvar = tf[i][j] + oneoverNMAX2*tempvar;
			/*
			tempvar = tf[i][j];
			error = discfunct(i,j,flag,tf,newfunct,tf[i][j]);
			while((fabs(error) > ROOT_ERROR || epsilon > ROOT_ERROR) && iterations < 50)
			{
				newvar = tempvar - 
					error/discdifunct(i,j,flag,tf,newfunct,tempvar);
				epsilon = fabs(newvar-tempvar);
				error = discfunct(i,j,flag,tf,newfunct,newvar);
				tempvar=newvar;
				iterations++;
			*/	/*if(iterations==99) printf("iterations at 49: (%d, %d) \n",i,j);*/
			/*}*/
			/* For Jacobi differencing */
			
			newfunct[i][j] = newvar;
			
			/* For Gauss-Seidel differencing */
			/*newfunct[i][j] = tf[i][j];
			tf[i][j] = newvar;*/
			/*if(xicoord(i)>0 && xicoord(i) <= 1.0- DELTA/2.0 &&
			  etacoord(j)>0 && etacoord(j) <= 1.0 - DELTA/2.0){
						cadj += exp(-2*newvar)*xicoord(i)*etacoord(j);
			}
			if(xicoord(i) <=1.0+DELTA/2.0 && xicoord(i) > 1.0-DELTA/2.0 && etacoord(j) <= 1.0-DELTA/2.){
				cadj += (-xicoord(i)+DELTA/2.0+1.0)/DELTA*exp(-2*newvar)*xicoord(i)*etacoord(j);
			}
			if(etacoord(j) <=1.0+DELTA/2 && etacoord(j) > 1.0-DELTA/2.0 && xicoord(i) <= 1.0-DELTA/2.){
				cadj += (-etacoord(j)+DELTA/2+1.0)/DELTA*exp(-2*newvar)*xicoord(i)*etacoord(j);
			}
			if(xicoord(i)<=1.0+DELTA/2.0 && xicoord(i) > 1.0-DELTA/2.0 &&
					etacoord(j) <=1.0+DELTA/2.0 && etacoord(j) > 1.0-DELTA/2.0 ){
				cadj += pow((-etacoord(j)+DELTA/2+1.0)/DELTA,2)*exp(-2*newvar)*xicoord(i)*etacoord(j);
			}*/

		}
	
	}
	
	
	/*	fshift = cadj * pow(DELTA,2);*/
	/*fshift = log(curvol/(cadj * pow(DELTA,2)))/2;*/

		tempvar = newfunct[NMAX/2][NMAX/2];

		for(i=1; i<=NMAX; i++){
			for(j=1; j<=NMAX; j++){
				/* For Jacobi differencing */
	
				temp = tf[i][j];			
				tf[i][j] = newfunct[i][j] - tempvar;
				newfunct[i][j] = temp;
				
				/* For Gauss-Seidel differencing */
				/*
				temp = tf[i][j] - fshift;
				tf[i][j] = temp;
				*/
/*				tf[i][j] = tf[i][j] + newfunct[NMAX][NMAX] -newvar;*/
			}
		}
		/**oldshift = (newfunct[NMAX][NMAX]-newvar)/oneoverNMAX2;*/
		return(fshift);
/*	}*/
}


double discfunct(int i, int j, int flag, double **tf)
{
	extern double DELTA;
	double DELTA4, test;
	DELTA4 = pow(DELTA,4);
	
if(flag == 1)
	return(
	 log(fabs((4*(tf[i+1][j]-2*tf[i][j]+tf[i-1][j])*(tf[i][j+1]-2*tf[i][j]+tf[i][j-1]) - 
	 1.0/16.0*pow(tf[i+1][j+1]-tf[i-1][j+1]-tf[i+1][j-1]+tf[i-1][j-1],2))/
	DELTA4))+ 2 * tf[i][j]
	);
else if(flag == 2){
	test = (2*(tf[i][-1 + j] - 2*tf[i][j] + tf[i][1 + j])*
        (tf[-1 + i][j] - 2*tf[i][j] + tf[1 + i][j] + 
          (-tf[-1 + i][j] + tf[1 + i][j])/(-1 + i)/2.) - 
       0.0625*pow(tf[-1 + i][-1 + j] - 
          tf[-1 + i][1 + j] - tf[1 + i][-1 + j] + 
          tf[1 + i][1 + j],2))/DELTA4;
	if(test <=0) printf("trouble at %d, %d \n", i,j);

	return( log(fabs(test)) + 2*tf[i][j]);
	}
else if(flag == 3)
	return(log(fabs((2*(tf[i][-1 + j] - 2*tf[i][j] + tf[i][1 + j] + 
          (-tf[i][-1 + j] + tf[i][1 + j])/(-1 + j)/2.)*
        (tf[-1 + i][j] - 2*tf[i][j] + tf[1 + i][j]) - 
       0.0625*pow(tf[-1 + i][-1 + j] - 
          tf[-1 + i][1 + j] - tf[1 + i][-1 + j] + 
          tf[1 + i][1 + j],2))/DELTA4)) + 2*tf[i][j]);
else{test = ((tf[i][-1 + j] - 2*tf[i][j] + tf[i][1 + j] + 
          (-tf[i][-1 + j] + tf[i][1 + j])/(-1 + j)/2.)*
        (tf[-1 + i][j] - 2*tf[i][j] + tf[1 + i][j] + 
          (-tf[-1 + i][j] + tf[1 + i][j])/(-1 + i)/2.) - 
       0.0625*pow(tf[-1 + i][-1 + j] - 
          tf[-1 + i][1 + j] - tf[1 + i][-1 + j] + 
          tf[1 + i][1 + j],2))/DELTA4;
	if(test<=0) printf("trouble at %d, %d \n", i,j);

	return( log(fabs(test)) + 2*tf[i][j]);
	};
}


void bcuint(y,y1,y2,y12,x1l,x1u,x2l,x2u,x1,x2,ansy,ansy1,ansy2)
double y[],y1[],y2[],y12[],x1l,x1u,x2l,x2u,x1,x2,*ansy,*ansy1,*ansy2;
{
	int i;
	double t,u,d1,d2,**c,**matrix();
	void bcucof(),nrerror(),free_matrix();

	c=matrix(1,4,1,4);
	d1=x1u-x1l;
	d2=x2u-x2l;
	bcucof(y,y1,y2,y12,d1,d2,c);
	if (x1u == x1l || x2u == x2l) {
		printf("badness \n");
		exit(1);
	};
	t=(x1-x1l)/d1;
	u=(x2-x2l)/d2;
	*ansy=(*ansy2)=(*ansy1)=0.0;
	for (i=4;i>=1;i--) {
		*ansy=t*(*ansy)+((c[i][4]*u+c[i][3])*u+c[i][2])*u+c[i][1];
		*ansy2=t*(*ansy2)+(3.0*c[i][4]*u+2.0*c[i][3])*u+c[i][2];
		*ansy1=u*(*ansy1)+(3.0*c[4][i]*t+2.0*c[3][i])*t+c[2][i];
	}
	*ansy1 /= d1;
	*ansy2 /= d2;
	free_matrix(c,1,4,1,4);
}


void bcucof(y,y1,y2,y12,d1,d2,c)
double y[],y1[],y2[],y12[],d1,d2,**c;
{
	static int wt[16][16]=
		{ 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
		0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,
		-3,0,0,3,0,0,0,0,-2,0,0,-1,0,0,0,0,
		2,0,0,-2,0,0,0,0,1,0,0,1,0,0,0,0,
		0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
		0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
		0,0,0,0,-3,0,0,3,0,0,0,0,-2,0,0,-1,
		0,0,0,0,2,0,0,-2,0,0,0,0,1,0,0,1,
		-3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,
		0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0,
		9,-9,9,-9,6,3,-3,-6,6,-6,-3,3,4,2,1,2,
		-6,6,-6,6,-4,-2,2,4,-3,3,3,-3,-2,-1,-1,-2,
		2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0,
		0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0,
		-6,6,-6,6,-3,-3,3,3,-4,4,2,-2,-2,-2,-1,-1,
		4,-4,4,-4,2,2,-2,-2,2,-2,-2,2,1,1,1,1};
	int l,k,j,i;
	double xx,d1d2,cl[16],x[16];

	d1d2=d1*d2;
	for (i=1;i<=4;i++) {
		x[i-1]=y[i];
		x[i+3]=y1[i]*d1;
		x[i+7]=y2[i]*d2;
		x[i+11]=y12[i]*d1d2;
	}
	for (i=0;i<=15;i++) {
		xx=0.0;
		for (k=0;k<=15;k++) xx += wt[i][k]*x[k];
		cl[i]=xx;
	}
	l=0;
	for (i=1;i<=4;i++)
		for (j=1;j<=4;j++) c[i][j]=cl[l++];
}

double **matrix(int imin,int imax,int jmin,int jmax){

int j;
double **tf;

tf = calloc(imax+1, sizeof(double *));
for(j=0; j<= jmax; j++){

 tf[j] = (double *) calloc(jmax+1, sizeof(double));
 if(tf[j]==NULL){
	printf("Memory allocation for grid failed. \n");
	exit(1);
 }
 }
 return(tf);
}

void free_matrix(double **c, int imin,int imax,int jmin,int jmax){

int j;
for(j=0; j<= jmax; j++){
 free(c[j]);
 }
 
 free(c);
}