/*	Copyright (c) 1997 Lucent Technologies	*/
/*	  All Rights Reserved 	*/

/*	THIS IS UNPUBLISHED PROPRIETARY SOURCE CODE OF LUCENT.	*/
/*	The copyright notice above does not evidence any   	*/
/*	actual or intended publication of such source code.	*/

/* sample funcadd */

#include "stdlib.h"	/* for atoi */
#include "math.h"
#include "funcadd.h"	/* includes "stdio1.h" */
#include "myalloc.h"

/* #define HUGE_VAL 1.0e+300 */

/* Sample function kth optionally prints its arguments using
 * a variant of stdio.h supplied by funcadd.h. */

#ifdef __cplusplus
/* The #ifdef __cplusplus lines allow this to be compiled
 * either by an ANSI C compiler or by a C++ compiler.
 */
extern "C" {
#endif

static void showvec( real *x, int n, AmplExports *ae);

static void showmat( real **A, int n, AmplExports *ae);

static void showmat2( real **A, int n, AmplExports *ae);

static void mx( real **A, real *x, real *y, int n);

static double dotprod( double *x, double *y, int n);

static int argmax( real *x, int n);

real erf (real x);
static real gammp (real a, real x);
static void gser(real *gamser, real a, real x, real *gln);
static void gcf( real *gammcf, real a, real x, real *gln);

 static real
ginv(arglist *al)	/* generalized inverse of a single argument */
{
	real x = al->ra[0];
	x = x ? 1./x : 0.;
	if (al->derivs) {
		*al->derivs = -x*x;
		if (al->hes)
			*al->hes = -2.*x * *al->derivs;
		}
	return x;
	}

 static char *
sginv(arglist *al)	/* character-valued version of ginv */
{
	static char buf[32];
	AmplExports *ae = al->AE;	/* for sprintf */
	real x = al->ra[0];
	sprintf(buf, "x%.g", x ? 1./x : 0);
	return buf;
	}

 static real
myhypot(arglist *al)	/* myhypot(x,y) = sqrt(x*x + y*y) */
{
	real *d, *h, rv, x, x0, y, y0;

	x = x0 = al->ra[0];
	y = y0 = al->ra[1];

	if (x < 0.)
		x = -x;
	if (y < 0.)
		y = -y;
	rv = x;
	if (x < y) {
		rv = y;
		y = x;
		x = rv;
		}
	if (rv) {
		y /= x;
		rv = x * sqrt(1. + y*y);
		if (d = (real *)al->derivs) {
			d[0] = x0 / rv;
			d[1] = y0 / rv;
			if (h = al->hes) {
				h[0] =  d[1]*d[1] / rv;
				h[1] = -d[0]*d[1] / rv;
				h[2] =  d[0]*d[0] / rv;
				}
			}
		}
	else if (d = (real *)al->derivs) {
		d[0] = d[1] = 0;
		if (h = al->hes)
			h[0] = h[1] = h[2] = 0;
		}
	return rv;
	}

 static real
ncall(void)		/* returns it's invocation count */
{ static real x; return ++x; }

 static real
mean(register arglist *al)	/* mean of arbitrarily many arguments */
{
	real x, z;
	real *d, *de, *ra;
	int *at, i, j, n, nreal;
	char *se, *sym;
	AmplExports *ae = al->AE;

	n = al->n;
	if (n <= 0) return 0;
	at = al->at;
	ra = al->ra;
	d = (real *)al->derivs;
	x = 0.0;
	nreal = 0;
	for(i=0; i<n; i++) {
		j = at[i];
		if (j >= 0) {
			x += ra[j];
			++nreal;
		} else {
			x += z = strtod(sym = al->sa[-(j+1)], &se);
			if (*se) {
				fprintf(Stderr,
				"mean treating arg %d = \"%s\" as %.g\n",
					i+1, sym, z);
				fflush(Stderr);
			}
		}
	}
	if (d) {
		z = 1.0/n;
		for (i=0; i<nreal; i++) {d[i] = z;}
		/* The Hessian is == 0 and, if needed, has been */
		/* initialized to 0. */
	}
	return x/n;
}

 static real
mineig(register arglist *al)	/* minimum eigenvalue of a symmetric matrix */
{
	real z, lambdak, normx, normy, shift=0, cosine, ydotEkk, 
		lambda0=HUGE_VAL, hm, num;
	real *e0, *ek, *d, *h, *ra;
	int *at, i, ii, j, jj, k, kk, m, n, narg, cnt=0, i0;
	char *se, *sym;
	AmplExports *ae = al->AE;
	static real **A=NULL, **E=NULL;
	static real *lambda=NULL, *x=NULL, *y=NULL, *y0=NULL;

	/*
	srand48(0);
	*/

	narg = al->n;
	if (narg <= 0) return 0;
	n = (int)(sqrt(1+8*narg)-1)/2;
	if (n*(n+1)/2 != narg) {
	     fprintf(Stderr, 
	     "should be mineig({i in 1..n, j in 1..n: j>=i} A[i,j])\n");
	     fflush(Stderr);
	}
	at = al->at;
	ra = al->ra;
	d = (real *)al->derivs;
	h = al->hes;

	/*
	if (A == NULL) {
	    CALLOC( A, n, real * );
	} else {
	    REALLOC( A, n, real * );
	}
	for (i=0; i<n; i++) REALLOC( A[i], n-i, real );
	*/
	if (A == NULL) {
	    REALLOC( A, n, real *);
	    A[0] = NULL;
	} else {
	    REALLOC( A, n, real *);
	}
	/*
	REALLOC( A[0], n*(n+1)/2, real);
	for (i=0, k=0; i<n; k+=n-i, i++) {
	    A[i] = A[0]+k;
	}
	*/
	REALLOC( A[0], n*n, real);
	for (i=0, k=0; i<n; k+=n, i++) A[i] = A[0]+k;

	if (E == NULL) {
	    REALLOC( E, n, real *);
	    E[0] = NULL;
	} else {
	    REALLOC( E, n, real *);
	}
	REALLOC( E[0], n*n, real);
	for (i=0, k=0; i<n; k+=n, i++) E[i] = E[0]+k;

	REALLOC( lambda, n, real );

	REALLOC( x, n, real );

	REALLOC( y, n, real );

	/*
	REALLOC( y0, n, real ); for (i=0; i<n; i++) y0[i] = drand48(); 
	*/
	REALLOC( y0, n, real ); for (i=0; i<n; i++) y0[i] = 0.5;

	k = 0;
	for(i=0; i<n; i++) {
	    for (j=i; j<n; j++) {
		jj = at[k];
fprintf(Stderr, "i,j,k,jj: %3d %3d %3d %3d\n",i,j,k,jj); fflush(Stderr);
		if (jj >= 0) {
			A[i][j] = ra[jj];
		} else {
			A[i][j] = z = strtod(sym = al->sa[-(jj+1)], &se);
			if (*se) {
				fprintf(Stderr,
				"mineig treating arg %d = \"%s\" as %.g\n",
					i+1, sym, z);
				fflush(Stderr);
			}
		}
		k++;
	    }
	}
	/* showmat(A,n,ae); */

	for (i=0; i<n; i++) y[i]  = y0[i];
	do {
	    normy = sqrt(dotprod(y,y,n));
	    for (i=0; i<n; i++) { y[i] /= normy; }
	    cosine = dotprod(x,y,n);
	    for (i=0; i<n; i++) { x[i] = y[i]; }
	    mx( A, x, y, n);
	    i = argmax(y,n);
	    lambda[0] = y[i]/x[i];
	    cnt++;
	} while ( fabs(cosine) < 0.99999999999 );

	/*
	fprintf(Stderr, "lambda = %10.4f \n", lambda[0]);
	if (lambda[0] >= 0) {shift = lambda[0]/2 + y[0];} else {shift = y[0];}
	*/
	shift = y0[0];
	fprintf(Stderr, "shift = %f \n", shift);

	for (i=0; i<n; i++) A[i][i] -= shift;

	for (k=0; k<n; k++) {
	    for (i=0; i<n; i++) y[i] = y0[i];
	    /*
	    for (kk=0; kk<k; kk++) {
		ydotEkk = dotprod(y, E[kk], n);
		for (i=0; i<n; i++) y[i] -= ydotEkk*E[kk][i];
	    }
	    */
	    do {
		normy = sqrt(dotprod(y,y,n));
	        for (i=0; i<n; i++) { y[i] /= normy; }
	        /* showvec(y,n,ae); */
		cosine = dotprod(x,y,n);
	        for (i=0; i<n; i++) { x[i] = y[i]; }
		mx( A, x, y, n);
	        /* showvec(y,n,ae); */
		i = argmax(y,n);
	        lambdak = y[i]/x[i];
		cnt++;
	    } while ( fabs(cosine) < 0.9999999999 );

	    lambda[k] = lambdak + shift;
	    for (i=0; i<n; i++) E[k][i] = x[i];

	    if (lambda[k] < lambda0) {
		i0 = k;
		lambda0 = lambda[k];
	    }

	    fprintf(Stderr, "%9.4f %9.4f: ", lambda[k], lambdak);
	    for (i=0; i<n; i++) { fprintf(Stderr, "%7.4f ", E[k][i]); }
	    fprintf(Stderr, "\n");

	    for (i=0; i<n; i++) {
		for (j=i; j<n; j++) {
		    A[i][j] -= lambdak*x[i]*x[j];
		}
	    }

	    /* showmat(A,n,ae); */
	}
	fprintf(Stderr, "\n");


	for(i=0; i<n; i++) {
	    for (j=i; j<n; j++) {
		A[i][j] = 0;
		for (k=0; k<n; k++) {
		    A[i][j] += lambda[k]*E[k][i]*E[k][j];
		}
	    }
	}
	showmat(A,n,ae);

	/*
	for (i=0; i<n; i++) {
	    for (j=i; j<n; j++) {
		fprintf(Stderr, "%2d %2d %10.4f \n", i, j,
		dotprod(E[i],E[j],n));
	    }
	}
	*/

	lambda0 = lambda[i0];
	e0 = E[i0];
	if (d) {
		m = 0;
		for (i=0; i<n; i++) {
		    for (j=i; j<n; j++) {
		        d[m] = 2*e0[i]*e0[j];
			if (i==j) d[m] /= 2;
			m++;
		    }
		}
	}
	if (h) {
	    m = 0;
	    for (i=0; i<n; i++) {
		for (j=i; j<n; j++) {
		    for (ii=i; ii<n; ii++) {
			for (jj=(ii==i)?j:ii; jj<n; jj++) {
			    hm = 0;
			    for (k=0; k<n; k++) {
				ek = E[k];
				if (k != i0) {
				    num = (e0[ i]*ek[ j] + e0[ j]*ek[ i])*
				          (e0[ii]*ek[jj] + e0[jj]*ek[ii]);
				    if ( i== j) num /= 2;
				    if (ii==jj) num /= 2;
				    hm += num/(lambda0-lambda[k]);
				}
			    }
			    h[m] = 2*hm;
			    m++;
			}
		    }
		}
	    }
	}
	return lambda0;
}

 static real
linear(arglist *al)	/* generalized inverse of a single argument */
{
	real x = al->ra[0];
	if (al->derivs) {
		*al->derivs = 1;
		if (al->hes)
			*al->hes = 0;
		}
	return x;
}

static real R_real, R_imag;
static real lx, ly, lz, qx, qy, qz, zmin, zmax;

static void
R(arglist *al) /* inverse of 3D spectral density on a spherical shell */
{
	int j, n=1000; 
/*	int j, n=20; */

	real pi = 4*atan(1);
	real lx, ly, lz, qx, qy, qz, zmin, zmax, z;
	real z_perp, 
	     lq, 	     /* component of l in direction of q */
	     lox, loy, loz,  /* component of l orthogonal to q */
	     l_mag, q_mag;

	lx = al->ra[0];
	ly = al->ra[1];
	lz = al->ra[2];
	qx = al->ra[3];
	qy = al->ra[4];
	qz = al->ra[5];
	zmin = al->ra[6];
	zmax = al->ra[7];

	q_mag = sqrt(qx*qx + qy*qy + qz*qz);
	qx /= q_mag;
	qy /= q_mag;
	qz /= q_mag;

	R_real=0;
	R_imag=0;

	for (j=0; j<n; j++) {
	    z = zmin + (zmax-zmin)*(j+0.5)/n;
	    z_perp = sqrt(1-z*z);
	    lq = lx*qx + ly*qy + lz*qz;
	    lox = lx - lq*qx;
	    loy = ly - lq*qy;
	    loz = lz - lq*qz;
	    l_mag = sqrt( lox*lox + loy*loy + loz*loz );
	    R_real += j0(2*pi*z_perp*l_mag)*cos(-2*pi*z*lq);
	    R_imag += j0(2*pi*z_perp*l_mag)*sin(-2*pi*z*lq);
	}
	R_real *= 2*pi*(zmax-zmin)/n;
	R_imag *= 2*pi*(zmax-zmin)/n;

	/*
	if (al->derivs) {
		*al->derivs = 1;
		if (al->hes)
			*al->hes = 0;
		}
	*/
}

 static real
C(arglist *al) /* inverse of 3D spectral density on a spherical shell */
{
	if (lx != al->ra[0] ||
	    ly != al->ra[1] ||
	    lz != al->ra[2] ||
	    qx != al->ra[3] ||
	    qy != al->ra[4] ||
	    qz != al->ra[5] ||
	    zmin != al->ra[6] ||
	    zmax != al->ra[7] ) {
		R(al);
	    }

	return R_real;
}

 static real
S(arglist *al) /* inverse of 3D spectral density on a spherical shell */
{
	if (lx != al->ra[0] ||
	    ly != al->ra[1] ||
	    lz != al->ra[2] ||
	    qx != al->ra[3] ||
	    qy != al->ra[4] ||
	    qz != al->ra[5] ||
	    zmin != al->ra[6] ||
	    zmax != al->ra[7] ) {
		R(al);
	    }

	return R_imag;
}

 static real
kth_eig(register arglist *al)	/* minimum kth of a symmetric matrix */
{
	real z, lambdak, normx, normy, shift=0, cosine, ydotEkk, 
		lambda0=HUGE_VAL, hm, num;
	real *e0, *ek, *d, *h, *ra;
	int *at, i, ii, j, jj, k, kk, m, n, N, cnt=0, i0;
	char *se, *sym;
	AmplExports *ae = al->AE;
	static real **A=NULL, **E=NULL;
	static real *lambda=NULL, *x=NULL, *y=NULL, *y0=NULL;

	/*
	srand48(0);
	*/

	N = al->n - 1;
	if (N <= 0) return 0;
	n = (int)(sqrt(1+8*N)-1)/2;
	if (n*(n+1)/2 != N) {
	     fprintf(Stderr, 
	     "should be kth_eig(k, {i in 1..n, j in 1..n: j>=i} A[i,j])\n");
	     fflush(Stderr);
	}
	at = al->at;
	ra = al->ra;
	d = (real *)al->derivs;
	h = al->hes;

	i0 = (int)al->ra[0];

	/*
	if (A == NULL) {
	    CALLOC( A, n, real * );
	} else {
	    REALLOC( A, n, real * );
	}
	for (i=0; i<n; i++) REALLOC( A[i], n-i, real );
	*/
	if (A == NULL) {
	    REALLOC( A, n, real *);
	    A[0] = NULL;
	} else {
	    REALLOC( A, n, real *);
	}
	/*
	REALLOC( A[0], n*(n+1)/2, real);
	for (i=0, k=0; i<n; k+=n-i, i++) {
	    A[i] = A[0]+k;
	}
	*/
	REALLOC( A[0], n*n, real);
	for (i=0, k=0; i<n; k+=n, i++) A[i] = A[0]+k;

	if (E == NULL) {
	    REALLOC( E, n, real *);
	    E[0] = NULL;
	} else {
	    REALLOC( E, n, real *);
	}
	REALLOC( E[0], n*n, real);
	for (i=0, k=0; i<n; k+=n, i++) E[i] = E[0]+k;

	REALLOC( lambda, n, real );

	REALLOC( x, n, real );

	REALLOC( y, n, real );

	if (y0 == NULL) {
	    /*
	    REALLOC( y0, n, real ); for (i=0; i<n; i++) y0[i] = drand48(); 
	    */
	    REALLOC( y0, n, real ); for (i=0; i<n; i++) y0[i] = 0.5; 
	} else {
	    REALLOC( y0, n, real ); 
	}

	k = 1;
	for(i=0; i<n; i++) {
	    for (j=i; j<n; j++) {
		jj = at[k];
		if (jj >= 0) {
			A[i][j] = ra[jj];
		} else {
			A[i][j] = z = strtod(sym = al->sa[-(jj+1)], &se);
			if (*se) {
				fprintf(Stderr,
				"kth_eig treating arg %d = \"%s\" as %.g\n",
					i+1, sym, z);
				fflush(Stderr);
			}
		}
		k++;
	    }
	}
	/* showmat(A,n,ae); */

	for (i=0; i<n; i++) y[i]  = y0[i];
	do {
	    normy = sqrt(dotprod(y,y,n));
	    for (i=0; i<n; i++) { y[i] /= normy; }
	    cosine = dotprod(x,y,n);
	    for (i=0; i<n; i++) { x[i] = y[i]; }
	    mx( A, x, y, n);
	    i = argmax(y,n);
	    lambda[0] = y[i]/x[i];
	    cnt++;
	} while ( fabs(cosine) < 0.99999999999 );

	shift = y0[0];

	for (i=0; i<n; i++) A[i][i] -= shift;

	for (k=0; k<n; k++) {
	    for (i=0; i<n; i++) y[i] = y0[i];
	    do {
		normy = sqrt(dotprod(y,y,n));
	        for (i=0; i<n; i++) { y[i] /= normy; }
	        /* showvec(y,n,ae); */
		cosine = dotprod(x,y,n);
	        for (i=0; i<n; i++) { x[i] = y[i]; }
		mx( A, x, y, n);
	        /* showvec(y,n,ae); */
		i = argmax(y,n);
	        lambdak = y[i]/x[i];
		cnt++;
	    } while ( fabs(cosine) < 0.9999999999 );

	    lambda[k] = lambdak + shift;
	    for (i=0; i<n; i++) E[k][i] = x[i];

	    for (i=0; i<n; i++) {
		for (j=i; j<n; j++) {
		    A[i][j] -= lambdak*x[i]*x[j];
		}
	    }

	    /* showmat(A,n,ae); */
	}

	for(i=0; i<n; i++) {
	    for (j=i; j<n; j++) {
		A[i][j] = 0;
		for (k=0; k<n; k++) {
		    A[i][j] += lambda[k]*E[k][i]*E[k][j];
		}
	    }
	}
	/* showmat(A,n,ae); */

	i0--;	/* convert from 1..n to 0..n-1 */
	lambda0 = lambda[i0];
	e0 = E[i0];
	if (d) {
		m = 0;
		for (i=0; i<n; i++) {
		    for (j=i; j<n; j++) {
		        d[m] = 2*e0[i]*e0[j];
			if (i==j) d[m] /= 2;
			m++;
		    }
		}
	}
	if (h) {
	    m = 0;
	    for (i=0; i<n; i++) {
		for (j=i; j<n; j++) {
		    for (ii=i; ii<n; ii++) {
			for (jj=(ii==i)?j:ii; jj<n; jj++) {
			    hm = 0;
			    for (k=0; k<n; k++) {
				ek = E[k];
				if (k != i0) {
				    num = (e0[ i]*ek[ j] + e0[ j]*ek[ i])*
				          (e0[ii]*ek[jj] + e0[jj]*ek[ii]);
				    if ( i== j) num /= 2;
				    if (ii==jj) num /= 2;
				    hm += num/(lambda0-lambda[k]);
				}
			    }
			    h[m] = 2*hm;
			    m++;
			}
		    }
		}
	    }
	}
	return lambda0;
}

 static real
kth_diag(register arglist *al)	/* kth diagonal element of D in LDL^T */
				/* must be called in order: 1,2,...,n */
{
	real z, diagval;
	real *d, *h, *ra;
	int *at, i, ii, j, jj, k, m, n, N;
	static int k0=-1;
	char *se, *sym;
	AmplExports *ae = al->AE;
	static real **A=NULL, Ident, **A0=NULL;
	static int convex_flag=1;

	N = al->n;
	if (N <= 2) return 0;
	n  = (int)al->ra[N-1];
	k0 = (int)al->ra[N-2];
	k0--;

	if (n*(n+1)/2 != N-2) {
	     fprintf(Stderr, 
	     "should be kth_diag({i in 1..n, j in 1..n: j>=i} A[i,j],k,n)\n");
	     fflush(Stderr);
	}
	at = al->at;
	ra = al->ra;
	d = (real *)al->derivs;
	h = al->hes;

/*
fprintf(Stderr,"IN KTH_DIAG \n");
fprintf(Stderr,"k0 = %3d \n", k0);
*/

	if (A == NULL) {
	    REALLOC( A, n, real *);
	    A[0] = NULL;
	} else {
	    REALLOC( A, n, real *);
	}
	REALLOC( A[0], 2*n*n, real);
	for (i=0, k=0; i<n; k+=2*n, i++) A[i] = A[0]+k;

	if (A0 == NULL) {
	    REALLOC( A0, n, real *);
	    A0[0] = NULL;
	} else {
	    REALLOC( A0, n, real *);
	}
	REALLOC( A0[0], n*n, real);
	for (i=0, k=0; i<n; k+=n, i++) A0[i] = A0[0]+k;

	if (k0 == 0) {
	    k = 0;
	    for(i=0; i<n; i++) {
	        for (j=i; j<n; j++) {
		    jj = at[k];
		    if (jj >= 0) {
			A[i][j] = ra[jj];
		    } else {
			A[i][j] = z = strtod(sym = al->sa[-(jj+1)], &se);
			if (*se) {
				fprintf(Stderr,
				"kth_diag treating arg %d = \"%s\" as %.g\n",
					i+1, sym, z);
				fflush(Stderr);
			}
		    }
		    A[j][i] = A[i][j];
		    k++;
	        }
	        for(j=n; j<2*n; j++) {
		    A[i][j] = 0;
	        }
		A[i][n+i] = 1;
	    }
	    for (i=0; i<n; i++) {
		for (j=0; j<n; j++) {
		    A0[i][j] = A[i][j];
		}
	    }
	    /*
	    showmat2(A,n,ae); 
	    */
	}

/*fprintf(Stderr,"Gradient \n");*/
	if (d) {
		m = 0;
		for (i=0; i<n; i++) {
		    for (j=i; j<n; j++) {
		        if (i<k0 && j<k0) {
			    d[m] = A[i][k0]*A[j][k0];
			    if (j!=i) d[m] *= 2;
		        } else 
		        if (i<k0 && j==k0) {
			    d[m] = -2*A[i][k0];
		        } else 
		        if (i==k0 && j==k0) {
			    d[m] = 1;
		        } else {
			    d[m] = 0;
		        }
/*
fprintf(Stderr,"%2d %2d %10f \n", i,j,d[m]);
*/
			m++;
		    }
		}
	}
/*
fprintf(Stderr,"Hessian \n");
*/
	if (h) {
	    m = 0;
	    for (ii=0; ii<n; ii++) {
		for (jj=ii; jj<n; jj++) {
		    for (i=0; i<=ii; i++) {
			for (j=i; j<=(i==ii?jj:n-1); j++) {
			    if (i<k0 && j<k0 && ii<k0 && jj<k0) {
			        if ( i!=j && ii!=jj) {
			            h[m] = - A[ii][k0]*A[jj][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+ii]*A[jj][k0]

			                   - A[ii][k0]*A[jj][n+j ]*A[i ][k0]
				           - A[j ][k0]*A[i ][n+ii]*A[jj][k0]

			                   - A[jj][k0]*A[ii][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+jj]*A[ii][k0]

			                   - A[jj][k0]*A[ii][n+j ]*A[i ][k0]
				           - A[j ][k0]*A[i ][n+jj]*A[ii][k0];
				} else 
				if ( i!=j && ii==jj) {
			            h[m] = - A[ii][k0]*A[jj][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+ii]*A[jj][k0]

			                   - A[ii][k0]*A[jj][n+j ]*A[i ][k0]
				           - A[j ][k0]*A[i ][n+ii]*A[jj][k0];
				} else 
				if ( i==j && ii!=jj) {
			            h[m] = - A[ii][k0]*A[jj][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+ii]*A[jj][k0]

			                   - A[jj][k0]*A[ii][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+jj]*A[ii][k0];
				} else 
				if ( i==j && ii==jj) {
			            h[m] = - A[ii][k0]*A[jj][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+ii]*A[jj][k0];
				}
			    } else 
			    if (i<k0 && j<k0 && ii<k0 && jj==k0) {
				h[m] = A[i][k0]*A[j][n+ii]
				     + A[j][k0]*A[ii][n+i];
				if (i!=j) h[m] *= 2;
			    } else 
			    if (i<k0 && j==k0 && ii<k0 && jj<k0) {
				h[m] = A[ii][k0]*A[jj][n+i]
				     + A[jj][k0]*A[i][n+ii];
				if (ii!=jj) h[m] *= 2;
			    } else 
			    if (i<k0 && j==k0 && ii<k0 && jj==k0) {
				h[m] = - A[ii][n+i ] - A[i ][n+ii];
			    } else {
				h[m] = 0;
			    }
/*
fprintf(Stderr,"%2d %2d %2d %2d %10f \n", i,j,ii,jj,h[m]);
*/
			    m++;
			}
		    }
		}
	    }
	}

	for (i=0; i<k0; i++) {
	    for (j=k0; j<n+k0; j++) {
		A[k0][j] -= A[k0][i]*A[i][j];
	    }
	    A[k0][i] = 0;
	}
	diagval = A[k0][k0];
/*
fprintf(Stderr,"val=%10f \n", diagval);
*/
	for (j=k0; j<=n+k0; j++) {
	    A[k0][j] /= diagval;
	}
	for (i=0; i<k0; i++) {
	    for (j=k0+1; j<=n+k0; j++) {
		A[i][j] -= A[i][k0]*A[k0][j];
	    }
	    A[i][k0] = 0;
	}

	/*
	fprintf(Stderr,"k = %3d \n", k0);
	showmat2(A,n,ae);
	fprintf(Stderr,"\n");
	fflush(Stderr);
	*/

	if (k0==n-1) {
	    for (i=0; i<n; i++) {
		for (j=0; j<n; j++) {
		    Ident = 0;
		    for (k=0; k<n; k++) {
			Ident += A0[i][k]*A[k][n+j];
		    }
		    if ( (i==j && fabs(Ident - 1) > 0.001)
		       ||(i!=j && fabs(Ident) > 0.001) ) {
			   fprintf(Stderr,"error in Gaussian elimination \n");
		       }
		}
	    }
	}

	if (diagval < -1.0e-12) {convex_flag = 0;}
	if (k0 == n-1) {
	    if (convex_flag == 0) {
		fprintf(Stderr,"matrix X not positive semidefinite \n");
	    }
	    convex_flag = 1;
	}

	return diagval;
}

 static real
kth_diag_eps(register arglist *al)	/* kth diagonal element of D in LDL^T */
				/* must be called in order: 1,2,...,n */
{
	real z, diagval;
	real *d, *h, *ra;
	int *at, i, ii, j, jj, k, m, n, N;
	static int k0=-1;
	char *se, *sym;
	AmplExports *ae = al->AE;
	static real **A=NULL, Ident, **A0=NULL;

	N = al->n;
	if (N <= 0) return 0;
	n = (int)(sqrt(1+8*N)-1)/2;
	if (n*(n+1)/2 != N) {
	     fprintf(Stderr, 
	     "should be kth_diag({i in 1..n, j in 1..n: j>=i} A[i,j])\n");
	     fflush(Stderr);
	}
	at = al->at;
	ra = al->ra;
	d = (real *)al->derivs;
	h = al->hes;

	k0++;  if (k0==n) k0 = 0;

/*
fprintf(Stderr,"IN KTH_DIAG \n");
fprintf(Stderr,"k0 = %3d \n", k0);
*/

	if (A == NULL) {
	    REALLOC( A, n, real *);
	    A[0] = NULL;
	} else {
	    REALLOC( A, n, real *);
	}
	REALLOC( A[0], 2*n*n, real);
	for (i=0, k=0; i<n; k+=2*n, i++) A[i] = A[0]+k;

	if (A0 == NULL) {
	    REALLOC( A0, n, real *);
	    A0[0] = NULL;
	} else {
	    REALLOC( A0, n, real *);
	}
	REALLOC( A0[0], n*n, real);
	for (i=0, k=0; i<n; k+=n, i++) A0[i] = A0[0]+k;

	if (k0 == 0) {
	    k = 0;
	    for(i=0; i<n; i++) {
	        for (j=i; j<n; j++) {
		    jj = at[k];
		    if (jj >= 0) {
			A[i][j] = ra[jj];
		    } else {
			A[i][j] = z = strtod(sym = al->sa[-(jj+1)], &se);
			if (*se) {
				fprintf(Stderr,
				"kth_diag treating arg %d = \"%s\" as %.g\n",
					i+1, sym, z);
				fflush(Stderr);
			}
		    }
		    A[j][i] = A[i][j];
		    k++;
	        }
	        for(j=n; j<2*n; j++) {
		    A[i][j] = 0;
	        }
		A[i][n+i] = 1;

		A[i][i] += 1.0e-6;
	    }
	    for (i=0; i<n; i++) {
		for (j=0; j<n; j++) {
		    A0[i][j] = A[i][j];
		}
	    }
	    /*
	    showmat2(A,n,ae); 
	    */
	}

/*
fprintf(Stderr,"Gradient \n");
*/
	if (d) {
		m = 0;
		for (i=0; i<n; i++) {
		    for (j=i; j<n; j++) {
		        if (i<k0 && j<k0) {
			    d[m] = A[i][k0]*A[j][k0];
			    if (j!=i) d[m] *= 2;
		        } else 
		        if (i<k0 && j==k0) {
			    d[m] = -2*A[i][k0];
		        } else 
		        if (i==k0 && j==k0) {
			    d[m] = 1;
		        } else {
			    d[m] = 0;
		        }
/*
fprintf(Stderr,"%2d %2d %10f \n", i,j,d[m]);
*/
			m++;
		    }
		}
	}
/*
fprintf(Stderr,"Hessian \n");
*/
	if (h) {
	    m = 0;
	    for (ii=0; ii<n; ii++) {
		for (jj=ii; jj<n; jj++) {
		    for (i=0; i<=ii; i++) {
			for (j=i; j<=(i==ii?jj:n-1); j++) {
			    if (i<k0 && j<k0 && ii<k0 && jj<k0) {
			        if ( i!=j && ii!=jj) {
			            h[m] = - A[ii][k0]*A[jj][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+ii]*A[jj][k0]

			                   - A[ii][k0]*A[jj][n+j ]*A[i ][k0]
				           - A[j ][k0]*A[i ][n+ii]*A[jj][k0]

			                   - A[jj][k0]*A[ii][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+jj]*A[ii][k0]

			                   - A[jj][k0]*A[ii][n+j ]*A[i ][k0]
				           - A[j ][k0]*A[i ][n+jj]*A[ii][k0];
				} else 
				if ( i!=j && ii==jj) {
			            h[m] = - A[ii][k0]*A[jj][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+ii]*A[jj][k0]

			                   - A[ii][k0]*A[jj][n+j ]*A[i ][k0]
				           - A[j ][k0]*A[i ][n+ii]*A[jj][k0];
				} else 
				if ( i==j && ii!=jj) {
			            h[m] = - A[ii][k0]*A[jj][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+ii]*A[jj][k0]

			                   - A[jj][k0]*A[ii][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+jj]*A[ii][k0];
				} else 
				if ( i==j && ii==jj) {
			            h[m] = - A[ii][k0]*A[jj][n+i ]*A[j ][k0]
				           - A[i ][k0]*A[j ][n+ii]*A[jj][k0];
				}
			    } else 
			    if (i<k0 && j<k0 && ii<k0 && jj==k0) {
				h[m] = A[i][k0]*A[j][n+ii]
				     + A[j][k0]*A[ii][n+i];
				if (i!=j) h[m] *= 2;
			    } else 
			    if (i<k0 && j==k0 && ii<k0 && jj<k0) {
				h[m] = A[ii][k0]*A[jj][n+i]
				     + A[jj][k0]*A[i][n+ii];
				if (ii!=jj) h[m] *= 2;
			    } else 
			    if (i<k0 && j==k0 && ii<k0 && jj==k0) {
				h[m] = - A[ii][n+i ] - A[i ][n+ii];
			    } else {
				h[m] = 0;
			    }
/*
fprintf(Stderr,"%2d %2d %2d %2d %10f \n", i,j,ii,jj,h[m]);
*/
			    m++;
			}
		    }
		}
	    }
	}

	for (i=0; i<k0; i++) {
	    for (j=k0; j<n+k0; j++) {
		A[k0][j] -= A[k0][i]*A[i][j];
	    }
	    A[k0][i] = 0;
	}
	diagval = A[k0][k0];
/*
fprintf(Stderr,"val=%10f \n", diagval);
*/
	for (j=k0; j<=n+k0; j++) {
	    A[k0][j] /= diagval;
	}
	for (i=0; i<k0; i++) {
	    for (j=k0+1; j<=n+k0; j++) {
		A[i][j] -= A[i][k0]*A[k0][j];
	    }
	    A[i][k0] = 0;
	}

	/*
	fprintf(Stderr,"k = %3d \n", k0);
	showmat2(A,n,ae);
	fprintf(Stderr,"\n");
	fflush(Stderr);
	*/

	if (k0==n-1) {
	    for (i=0; i<n; i++) {
		for (j=0; j<n; j++) {
		    Ident = 0;
		    for (k=0; k<n; k++) {
			Ident += A0[i][k]*A[k][n+j];
		    }
		    if ( (i==j && fabs(Ident - 1) > 0.001)
		       ||(i!=j && fabs(Ident) > 0.001) ) {
			   fprintf(Stderr,"error in Gaussian elimination \n");
		       }
		}
	    }
	}

	return diagval;
}

 static real
kth_diag2(register arglist *al)	/* kth diagonal element of D in LDL^T */
				/* must be called in order: 1,2,...,n */
{
	real z, diagval;
	real *d, *h, *ra;
	int *at, i, ii, j, jj, k, m, n, N;
	static int k0=-1;
	char *se, *sym;
	AmplExports *ae = al->AE;
	static real **A=NULL, Ident, **A0=NULL;

	N = al->n - 1;
	if (N <= 0) return 0;
	n = (int)(sqrt(N));
	if (n*n != N) {
	     fprintf(Stderr, 
	     "should be kth_diag2(k, {i in 1..n, j in 1..n} A[i,j])\n");
	     fflush(Stderr);
	}
	at = al->at;
	ra = al->ra;
	d = (real *)al->derivs;
	h = al->hes;

	k = (int)al->ra[0];
	k--;	/* convert from 1..n to 0..n-1 */
	if (k != k0+1 && k != 0) {
	     fprintf(Stderr, 
	     "k in kth_diag2 not in increasing order \n");
	     fflush(Stderr);
	}
	k0 = k;

	if (A == NULL) {
	    REALLOC( A, n, real *);
	    A[0] = NULL;
	} else {
	    REALLOC( A, n, real *);
	}
	REALLOC( A[0], 2*n*n, real);
	for (i=0, k=0; i<n; k+=2*n, i++) A[i] = A[0]+k;

	if (A0 == NULL) {
	    REALLOC( A0, n, real *);
	    A0[0] = NULL;
	} else {
	    REALLOC( A0, n, real *);
	}
	REALLOC( A0[0], n*n, real);
	for (i=0, k=0; i<n; k+=n, i++) A0[i] = A0[0]+k;

	if (k0 == 0) {
	    k = 1;
	    for(i=0; i<n; i++) {
	        for (j=0; j<n; j++) {
		    jj = at[k];
		    if (jj >= 0) {
			A[i][j] = ra[jj];
		    } else {
			A[i][j] = z = strtod(sym = al->sa[-(jj+1)], &se);
			if (*se) {
				fprintf(Stderr,
				"kth_diag2 treating arg %d = \"%s\" as %.g\n",
					i+1, sym, z);
				fflush(Stderr);
			}
		    }
		    k++;
	        }
	        for(j=n; j<2*n; j++) {
		    A[i][j] = 0;
	        }
		A[i][n+i] = 1;
	    }
	    for (i=0; i<n; i++) {
		for (j=0; j<n; j++) {
		    A0[i][j] = A[i][j];
		}
	    }
	    showmat2(A,n,ae);
	}

	if (d) {
		m = 0;
		for (i=0; i<n; i++) {
		    for (j=0; j<n; j++) {
		        if (i<k0 && j<k0) {
			    d[m] = A[i][k0]*A[j][k0];
		        } else 
		        if (i<k0 && j==k0) {
			    d[m] = -2*A[i][k0];
		        } else 
		        if (i==k0 && j<k0) {
			    d[m] = -2*A[j][k0];
		        } else 
		        if (i==k0 && j==k0) {
			    d[m] = 1;
		        } else {
			    d[m] = 0;
		        }
			m++;
		    }
		}
	}
	if (h) {
	    m = 0;
	    for (i=0; i<n; i++) {
		for (j=0; j<n; j++) {
		    for (ii=i; ii<n; ii++) {
			for (jj=(ii==i)?j:0; jj<n; jj++) {
			    if (i<k0 && j<k0 && ii<k0 && jj<k0) {
			        h[m] = - A[ii][k0]*A[jj][n+i ]*A[j ][k0]
				       - A[i ][k0]*A[j ][n+ii]*A[jj][k0];
			    } else 
			    if (i<k0 && j<k0 && ii<k0 && jj==k0) {
				h[m] = A[i][k0]*A[j][n+ii]
				     + A[j][k0]*A[ii][n+i];
			    } else 
			    if (i<k0 && j<k0 && ii==k0 && jj<k0) {
				h[m] = A[i][k0]*A[j][n+jj]
				     + A[j][k0]*A[jj][n+i];
			    } else 
			    if (i<k0 && j==k0 && ii<k0 && jj<k0) {
				h[m] = A[ii][k0]*A[jj][n+i]
				     + A[jj][k0]*A[i][n+ii];
			    } else 
			    if (i<k0 && j==k0 && ii<k0 && jj==k0) {
				h[m] = - A[ii][n+i ] - A[i ][n+ii];
			    } else 
			    if (i<k0 && j==k0 && ii==k0 && jj<k0) {
				h[m] = - A[jj][n+i ] - A[i ][n+jj];
			    } else 
			    if (i==k0 && j<k0 && ii<k0 && jj<k0) {
				h[m] = A[ii][k0]*A[jj][n+j]
				     + A[jj][k0]*A[j][n+ii];
			    } else 
			    if (i==k0 && j<k0 && ii<k0 && jj==k0) {
				h[m] = - A[ii][n+j ] - A[j ][n+ii];
			    } else 
			    if (i==k0 && j<k0 && ii==k0 && jj<k0) {
				h[m] = - A[jj][n+j ] - A[j ][n+jj];
			    } else {
				h[m] = 0;
			    }
			    m++;
			}
		    }
		}
	    }
	}

	for (i=0; i<k0; i++) {
	    for (j=k0; j<n+k0; j++) {
		A[k0][j] -= A[k0][i]*A[i][j];
	    }
	    A[k0][i] = 0;
	}
	diagval = A[k0][k0];
	for (j=k0; j<=n+k0; j++) {
	    A[k0][j] /= diagval;
	}
	for (i=0; i<k0; i++) {
	    for (j=k0+1; j<=n+k0; j++) {
		A[i][j] -= A[i][k0]*A[k0][j];
	    }
	    A[i][k0] = 0;
	}

	/*
	fprintf(Stderr,"k = %3d \n", k0);
	showmat2(A,n,ae);
	fprintf(Stderr,"\n");
	fflush(Stderr);
	*/

	if (k0==n-1) {
	    for (i=0; i<n; i++) {
		for (j=0; j<n; j++) {
		    Ident = 0;
		    for (k=0; k<n; k++) {
			Ident += A0[i][k]*A[k][n+j];
		    }
		    if ( (i==j && fabs(Ident - 1) > 0.00001)
		       ||(i!=j && fabs(Ident) > 0.00001) ) {
			   fprintf(Stderr,"error in Gaussian elimination \n");
		       }
		}
	    }
	}

	return diagval;
}

 static real
myerf(arglist *al)	/* 1/sqrt(2 pi) int_-infty^x e^{-t^2/2} dt */
{
	int j;
	double h = 1.0/128.0, x0 = -10;
	double h2 = h/2;
	double myerf=0;
	double pi = 4*atan(1);
	double rtp = sqrt(2*pi);
	double rt  = sqrt(2);
	double deriv;

	real x = al->ra[0];

	myerf = (erf(x/rt) + 1)/2;

	if (al->derivs) {
		deriv = exp(-x*x/2)/rtp;
		*al->derivs = deriv;
		if (al->hes)
			*al->hes = -x * deriv;
		}
	return myerf;
	}

 static real
stderf(arglist *al)	/* 2/sqrt(pi) int_0^x e^{-t^2} dt */
{
	int j;
	double h = 1.0/128.0, x0 = -10;
	double h2 = h/2;
	double stderf=0;
	double pi = 4*atan(1);
	double rp = sqrt(pi);
	double rt  = sqrt(2);
	double deriv;

	real x = al->ra[0];

	stderf = erf(x);

	if (al->derivs) {
		deriv = 2*exp(-x*x)/rp;
		*al->derivs = deriv;
		if (al->hes)
			*al->hes = -2*x * deriv;
		}
	return stderf;
	}

static void showvec( real *x, int n, AmplExports *ae)
{
	int i, j;

	for(i=0; i<n; i++) {
	    fprintf(Stderr, "%10.2e ", x[i]);
	}
	fprintf(Stderr, "\n");
	fflush(Stderr);
}

static void showmat( real **A, int n, AmplExports *ae)
{
	int i, j;

	for(i=0; i<n; i++) {
	    for (j=0; j<n; j++) {
		if (j>=i) {
		    fprintf(Stderr, "%10.4f ", A[i][j]);
		} else {
		    fprintf(Stderr, "           ");
		}
	    }
	    fprintf(Stderr, "\n");
	    fflush(Stderr);
	}
}

static void showmat2( real **A, int n, AmplExports *ae)
{
	int i, j;

	for(i=0; i<n; i++) {
	    for (j=0; j<2*n; j++) {
		    fprintf(Stderr, "%10.4f ", A[i][j]);
	    }
	    fprintf(Stderr, "\n");
	    fflush(Stderr);
	}
}

static void mx( real **A, real *x, real *y, int n)
{
	int i,j;

	for (i=0; i<n; i++) y[i] = 0;
	for (i=0; i<n; i++) {
	    y[i] += A[i][i]*x[i];
	    for (j=i+1; j<n; j++) {
		y[i] += A[i][j]*x[j];
		y[j] += A[i][j]*x[i];
	    }
	}
}

static real dotprod(         /* inner product between n-vectors x and y */
        real  *x,
        real  *y,
        int     n
)
{
        int i;
        real dotprod=0.0e0;

        for (i=0; i<n; i++) dotprod += x[i]*y[i];

        return dotprod;
}

static int argmax(         /* inner product between n-vectors x and y */
        real  *x,
        int     n
)
{
        int i, ii = -1;
        real max = -1.0e0;

        for (i=0; i<n; i++) if (fabs(x[i]) > max) {ii=i; max=fabs(x[i]);}

        return ii;
}

 static real
mymin(register arglist *al)	/* minimum eigenvalue of a symmetric matrix */
{
	real z, x0=HUGE_VAL;
	real *d, *h, *ra;
	int *at, i, j, jj, k, m, n, i0;
	char *se, *sym;
	AmplExports *ae = al->AE;
	static real *x=NULL;

	n = al->n;
	if (n <= 0) return 0;
	at = al->at;
	ra = al->ra;
	d = (real *)al->derivs;
	h = al->hes;

	REALLOC( x, n, real );

	for(i=0; i<n; i++) {
		jj = at[i];
		if (jj >= 0) {
			x[i] = ra[jj];
		} else {
			x[i] = z = strtod(sym = al->sa[-(jj+1)], &se);
			if (*se) {
				fprintf(Stderr,
				"mymin treating arg %d = \"%s\" as %.g\n",
					i+1, sym, z);
				fflush(Stderr);
			}
		}
	}

	for (k=0; k<n; k++) {
	    if (x[k] < x0) {
		i0 = k;
		x0 = x[k];
	    }
	}
/*
fprintf(Stderr, "min var = %d, val = %10.5f \n", i0, x0); fflush(Stderr);
*/

	if (d) {
		for (i=0; i<n; i++) { d[i] = 0; }
		d[i0] = 1;
	}
	if (h) {
	    m = 0;
	    for (i=0; i<n; i++) {
		for (j=i; j<n; j++) {
		    h[m] = 0;
		    m++;
		}
	    }
	}

	return x0;
}

 static real
nonneg(register arglist *al)	/* minimum eigenvalue of a symmetric matrix */
{
	real z, x0=HUGE_VAL;
	real *d, *h, *ra;
	int *at, i, j, jj, k, m, n, i0;
	char *se, *sym;
	AmplExports *ae = al->AE;
	static real *x=NULL;

	n = al->n - 1;
	if (n <= 0) return 0;
	at = al->at;
	ra = al->ra;
	d = (real *)al->derivs;
	h = al->hes;

	REALLOC( x, n+1, real );

	i0 = (int)al->ra[0];

	for(i=1; i<=n; i++) {
		jj = at[i];
		if (jj >= 0) {
			x[i] = ra[jj];
		} else {
			x[i] = z = strtod(sym = al->sa[-(jj+1)], &se);
			if (*se) {
				fprintf(Stderr,
				"nonneg treating arg %d = \"%s\" as %.g\n",
					i+1, sym, z);
				fflush(Stderr);
			}
		}
	}

	if (d) {
		for (i=0; i<n; i++) { d[i] = 0; }
		d[i0-1] = 1;
	}
	if (h) {
	    m = 0;
	    for (i=0; i<n; i++) {
		for (j=i; j<n; j++) {
		    h[m] = 0;
		    m++;
		}
	    }
	}

	return x[i0];
}

 static char *
kth(register arglist *al)	/* kth(k,a1,a2,...,an) return ak */
{
	int i, j, k, n;
	char *comma;
	static char buf[32];
	AmplExports *ae = al->AE;

	k = al->at[0] ? atoi(al->sa[0]) : (int)al->ra[0];
	n = al->n;
	if (k < 0) {
		fprintf(Stderr, "kth(");
		comma = "";
		for(i = 0; i < n; i++, comma = ", ")
			if ((j = al->at[i]) >= 0)
				fprintf(Stderr, "%s%.g", comma, al->ra[j]);
			else
				fprintf(Stderr, "%s%s", comma, al->sa[-(j+1)]);
		fprintf(Stderr, ")\n");
		fflush(Stderr);
		k = -k;
		}
	if (n <= 1 || k <= 0 || k >= n)
		return "";
	if ((j = al->at[k]) >= 0) {
		sprintf(buf, "%.g", al->ra[j]);
		return buf;
		}
	return al->sa[-(j+1)];
	}

 static real
junk(register arglist *al)	/* kth diagonal element of D in LDL^T */
				/* must be called in order: 1,2,...,n */
{
	real z, diagval;
	real *d, *h, *ra;
	int *at, i, j, jj, k, m, n, N;
	static int k0=-1;
	char *se, *sym;
	AmplExports *ae = al->AE;
	static real **A=NULL;

	N = al->n;
	if (N <= 0) return 0;
	n = (int)(sqrt(1+8*N)-1)/2;
	if (n*(n+1)/2 != N) {
	     fprintf(Stderr, 
	     "should be kth_diag({i in 1..n, j in 1..n: j>=i} A[i,j])\n");
	     fflush(Stderr);
	}
	at = al->at;
	ra = al->ra;
	d = (real *)al->derivs;
	h = al->hes;

	k0++;  if (k0==n) k0 = 0;

/*
fprintf(Stderr,"IN KTH_DIAG \n");
fprintf(Stderr,"k0 = %3d \n", k0);
*/

	if (A == NULL) {
	    REALLOC( A, n, real *);
	    A[0] = NULL;
	} else {
	    REALLOC( A, n, real *);
	}
	REALLOC( A[0], 2*n*n, real);
	for (i=0, k=0; i<n; k+=2*n, i++) A[i] = A[0]+k;

	if (k0 == 0) {
	    k = 0;
	    for(i=0; i<n; i++) {
	        for (j=i; j<n; j++) {
		    jj = at[k];
		    if (jj >= 0) {
			A[i][j] = ra[jj];
		    }
		    A[j][i] = A[i][j];
		    k++;
	        }
	        for(j=n; j<2*n; j++) {
		    A[i][j] = 0;
	        }
		A[i][n+i] = 1;
	    }
	    /**/
	    showmat2(A,n,ae); 
	    /**/
	}

	return 1;
}

static void alatalon(AmplExports *ae);
static double **aLat=NULL;
static double **aLon=NULL;
static int firstpass=1;

#define SS 0.5  /* degrees of lat/lon between data points */
#define leftlon   (-120)
#define rightlon  (-62.5)
#define bottomlat ( 20.5)
#define toplat    ( 50)

/* latitude component of wind velocity in deg/hr */
 static real
aLatf(register arglist *al)	
{
        int i, j;
	int m=1+(29.5/SS), n=1+(57.5/SS);
	real *d, *h, rv=0, lat, lat0, lon, lon0;
	AmplExports *ae = al->AE;
	real pi = 4*atan(1);

	if (firstpass) {
	    firstpass=0;
	    alatalon(ae);
	}

	lat = lat0 =  al->ra[0];
	lon = lon0 =  al->ra[1];  

	if (lat > 90) {lat = 180-lat;}
	if (lat < -90) {lat = -180+lat;}
	if (lon > 180) {lon -= 360;}
	if (lon < -180) {lon += 360;}
	i = floor((lat- bottomlat)/SS); if (i<1) {i=1;} if (i>=m-2) {i=m-3;}
	j = floor((lon- leftlon)/SS);   if (j<1) {j=1;} if (j>=n-2) {j=n-3;}

	rv = (lat-(SS*(i-1)+bottomlat))    *(lon-(leftlon+SS*(j-1)))     *aLat[i+1][j+1] + 
	     ((SS*(i+2)+bottomlat)-lat)    *(lon-(leftlon+SS*(j-1)))     *aLat[i  ][j+1] + 
	     (lat-(SS*(i-1)+bottomlat))    *((leftlon+SS*(j+2))-lon)     *aLat[i+1][j  ] + 
	     ((SS*(i+2)+bottomlat)-lat)    *((leftlon+SS*(j+2))-lon)     *aLat[i  ][j  ] +
	     (lat-(SS*i+bottomlat))        *((leftlon+SS*(j+1))-lon)     *aLat[i+2][j-1] + 
	     ((SS*(i+1)+bottomlat)-lat)    *((leftlon+SS*(j+1))-lon)     *aLat[i-1][j-1] +
	     ((SS*(i+1)+bottomlat)-lat)    *(lon-(leftlon+SS*j))         *aLat[i-1][j+2] + 
	     (lat-(SS*i+bottomlat))        *(lon-(leftlon+SS*j))         *aLat[i+2][j+2] + 
	     (lat-(SS*i+bottomlat))        *((leftlon+SS*(j+2))-lon)     *aLat[i+1][j-1] + 
	     ((SS*(i+1)+bottomlat)-lat)    *((leftlon+SS*(j+2))-lon)     *aLat[i  ][j-1] +
	     ((SS*(i+2)+bottomlat)-lat)    *((leftlon+SS*(j+1))-lon)     *aLat[i-1][j  ] +
	     ((SS*(i+2)+bottomlat)-lat)    *(lon-(leftlon+SS*j))         *aLat[i-1][j+1] + 
	     ((SS*(i+1)+bottomlat)-lat)    *(lon-(leftlon+SS*(j-1)))     *aLat[i  ][j+2] + 
	     (lat-(SS*i+bottomlat))        *(lon-(leftlon+SS*(j-1)))     *aLat[i+1][j+2] + 
	     (lat-(SS*(i-1)+bottomlat))    *(lon-(leftlon+SS*j))         *aLat[i+2][j+1] + 
	     (lat-(SS*(i-1)+bottomlat))    *((leftlon+SS*(j+1))-lon)     *aLat[i+2][j  ];

	rv /= 16*SS*SS;

	if (d = (real *)al->derivs) {
	       d[0] = (1      )    *(lon-(leftlon+SS*(j-1)))     *aLat[i+1][j+1] + 
	              (     -1)    *(lon-(leftlon+SS*(j-1)))     *aLat[i  ][j+1] + 
	              (1      )    *((leftlon+SS*(j+2))-lon)     *aLat[i+1][j  ] + 
	              (     -1)    *((leftlon+SS*(j+2))-lon)     *aLat[i  ][j  ] +
	              (1      )    *((leftlon+SS*(j+1))-lon)     *aLat[i+2][j-1] + 
	              (     -1)    *((leftlon+SS*(j+1))-lon)     *aLat[i-1][j-1] +
	              (     -1)    *(lon-(leftlon+SS*j))         *aLat[i-1][j+2] + 
	              (1      )    *(lon-(leftlon+SS*j))         *aLat[i+2][j+2] + 
	              (1      )    *((leftlon+SS*(j+2))-lon)     *aLat[i+1][j-1] + 
	              (     -1)    *((leftlon+SS*(j+2))-lon)     *aLat[i  ][j-1] +
	              (     -1)    *((leftlon+SS*(j+1))-lon)     *aLat[i-1][j  ] +
	              (     -1)    *(lon-(leftlon+SS*j))         *aLat[i-1][j+1] + 
	              (     -1)    *(lon-(leftlon+SS*(j-1)))     *aLat[i  ][j+2] + 
	              (1      )    *(lon-(leftlon+SS*(j-1)))     *aLat[i+1][j+2] + 
	              (1      )    *(lon-(leftlon+SS*j))         *aLat[i+2][j+1] + 
	              (1      )    *((leftlon+SS*(j+1))-lon)     *aLat[i+2][j  ];

	      d[0] /= 16*SS*SS;

	      d[1] = (lat-(SS*(i-1)+bottomlat))    *(1      )     *aLat[i+1][j+1] + 
	             ((SS*(i+2)+bottomlat)-lat)    *(1      )     *aLat[i  ][j+1] + 
	             (lat-(SS*(i-1)+bottomlat))    *(     -1)     *aLat[i+1][j  ] + 
	             ((SS*(i+2)+bottomlat)-lat)    *(     -1)     *aLat[i  ][j  ] +
	             (lat-(SS*i+bottomlat))        *(     -1)     *aLat[i+2][j-1] + 
	             ((SS*(i+1)+bottomlat)-lat)    *(     -1)     *aLat[i-1][j-1] +
	             ((SS*(i+1)+bottomlat)-lat)    *(1      )     *aLat[i-1][j+2] + 
	             (lat-(SS*i+bottomlat))        *(1      )     *aLat[i+2][j+2] + 
       	             (lat-(SS*i+bottomlat))        *(     -1)     *aLat[i+1][j-1] + 
	             ((SS*(i+1)+bottomlat)-lat)    *(     -1)     *aLat[i  ][j-1] +
	             ((SS*(i+2)+bottomlat)-lat)    *(     -1)     *aLat[i-1][j  ] +
	             ((SS*(i+2)+bottomlat)-lat)    *(1      )     *aLat[i-1][j+1] + 
	             ((SS*(i+1)+bottomlat)-lat)    *(1      )     *aLat[i  ][j+2] + 
	             (lat-(SS*i+bottomlat))        *(1      )     *aLat[i+1][j+2] + 
	             (lat-(SS*(i-1)+bottomlat))    *(1      )     *aLat[i+2][j+1] + 
	             (lat-(SS*(i-1)+bottomlat))    *(     -1)     *aLat[i+2][j  ];

		d[1] /= 16*SS*SS;

		if (h = al->hes) {
			h[0] =  0;
			h[1] = aLat[i+1][j+1]
			      -aLat[i  ][j+1]
			      -aLat[i+1][j  ]
			      +aLat[i  ][j  ]
			      -aLat[i+2][j-1]
			      +aLat[i-1][j-1]
			      -aLat[i-1][j+2]
			      +aLat[i+2][j+2]
			      -aLat[i+1][j-1]
			      +aLat[i  ][j-1]
			      +aLat[i-1][j  ]
			      -aLat[i-1][j+1]
			      -aLat[i  ][j+2]
			      +aLat[i+1][j+2]
			      +aLat[i+2][j+1]
			      -aLat[i+2][j  ];

			h[1] /= 16*SS*SS;
			h[2] =  0;
			}
		}

	return (rv);
	}

/* longitude component of wind velocity in deg/hr */
 static real
aLonf(register arglist *al)	
{
        int i, j;
	int m=1+(29.5/SS), n=1+(57.5/SS);
	real *d, *h, rv=0, lat, lon;
	real pi = 4*atan(1);
	AmplExports *ae = al->AE;

	if (firstpass) {
	    firstpass=0;
	    alatalon(ae);
	    }
/*
	    REALLOC(aLon, m, real *);
	    for (i=0; i<m; i++) {
		REALLOC(aLon[i], n, real);
		for (j=0; j<n; j++) {
		    aLon[i][j] = i-8*m/12;
		}
	    }
*/

	lat =  al->ra[0];
	lon =  al->ra[1];  

	if (lat > 90) {lat = 180-lat;}
	if (lat < -90) {lat = -180+lat;}
	if (lon > 180) {lon -= 360;}
	if (lon < -180) {lon += 360;}
	i = floor((lat-bottomlat)/SS); if (i<1) {i=1;} if (i>=m-2) {i=m-3;}
	j = floor((lon-leftlon)/SS);   if (j<1) {j=1;} if (j>=n-2) {j=n-3;}

	rv = (lat-(SS*(i-1)+bottomlat))    *(lon-(leftlon+SS*(j-1)))     *aLon[i+1][j+1] + 
	     ((SS*(i+2)+bottomlat)-lat)    *(lon-(leftlon+SS*(j-1)))     *aLon[i  ][j+1] + 
	     (lat-(SS*(i-1)+bottomlat))    *((leftlon+SS*(j+2))-lon)     *aLon[i+1][j  ] + 
	     ((SS*(i+2)+bottomlat)-lat)    *((leftlon+SS*(j+2))-lon)     *aLon[i  ][j  ] +
	     (lat-(SS*i+bottomlat))        *((leftlon+SS*(j+1))-lon)     *aLon[i+2][j-1] + 
	     ((SS*(i+1)+bottomlat)-lat)    *((leftlon+SS*(j+1))-lon)     *aLon[i-1][j-1] +
	     ((SS*(i+1)+bottomlat)-lat)    *(lon-(leftlon+SS*j))         *aLon[i-1][j+2] + 
	     (lat-(SS*i+bottomlat))        *(lon-(leftlon+SS*j))         *aLon[i+2][j+2] + 
	     (lat-(SS*i+bottomlat))        *((leftlon+SS*(j+2))-lon)     *aLon[i+1][j-1] + 
	     ((SS*(i+1)+bottomlat)-lat)    *((leftlon+SS*(j+2))-lon)     *aLon[i  ][j-1] +
	     ((SS*(i+2)+bottomlat)-lat)    *((leftlon+SS*(j+1))-lon)     *aLon[i-1][j  ] +
	     ((SS*(i+2)+bottomlat)-lat)    *(lon-(leftlon+SS*j))         *aLon[i-1][j+1] + 
	     ((SS*(i+1)+bottomlat)-lat)    *(lon-(leftlon+SS*(j-1)))     *aLon[i  ][j+2] + 
	     (lat-(SS*i+bottomlat))        *(lon-(leftlon+SS*(j-1)))     *aLon[i+1][j+2] + 
	     (lat-(SS*(i-1)+bottomlat))    *(lon-(leftlon+SS*j))         *aLon[i+2][j+1] + 
	     (lat-(SS*(i-1)+bottomlat))    *((leftlon+SS*(j+1))-lon)     *aLon[i+2][j  ];

	rv /= 16*SS*SS;

	if (d = (real *)al->derivs) {

	       d[0] = (1      )    *(lon-(leftlon+SS*(j-1)))     *aLon[i+1][j+1] + 
	              (     -1)    *(lon-(leftlon+SS*(j-1)))     *aLon[i  ][j+1] + 
	              (1      )    *((leftlon+SS*(j+2))-lon)     *aLon[i+1][j  ] + 
	              (     -1)    *((leftlon+SS*(j+2))-lon)     *aLon[i  ][j  ] +
	              (1      )    *((leftlon+SS*(j+1))-lon)     *aLon[i+2][j-1] + 
	              (     -1)    *((leftlon+SS*(j+1))-lon)     *aLon[i-1][j-1] +
	              (     -1)    *(lon-(leftlon+SS*j))         *aLon[i-1][j+2] + 
	              (1      )    *(lon-(leftlon+SS*j))         *aLon[i+2][j+2] + 
	              (1      )    *((leftlon+SS*(j+2))-lon)     *aLon[i+1][j-1] + 
	              (     -1)    *((leftlon+SS*(j+2))-lon)     *aLon[i  ][j-1] +
	              (     -1)    *((leftlon+SS*(j+1))-lon)     *aLon[i-1][j  ] +
	              (     -1)    *(lon-(leftlon+SS*j))         *aLon[i-1][j+1] + 
	              (     -1)    *(lon-(leftlon+SS*(j-1)))     *aLon[i  ][j+2] + 
	              (1      )    *(lon-(leftlon+SS*(j-1)))     *aLon[i+1][j+2] + 
	              (1      )    *(lon-(leftlon+SS*j))         *aLon[i+2][j+1] + 
	              (1      )    *((leftlon+SS*(j+1))-lon)     *aLon[i+2][j  ];

		d[0] /= 16*SS*SS;

	      d[1] = (lat-(SS*(i-1)+bottomlat))    *(1      )     *aLon[i+1][j+1] + 
	             ((SS*(i+2)+bottomlat)-lat)    *(1      )     *aLon[i  ][j+1] + 
	             (lat-(SS*(i-1)+bottomlat))    *(     -1)     *aLon[i+1][j  ] + 
	             ((SS*(i+2)+bottomlat)-lat)    *(     -1)     *aLon[i  ][j  ] +
	             (lat-(SS*i+bottomlat))        *(     -1)     *aLon[i+2][j-1] + 
	             ((SS*(i+1)+bottomlat)-lat)    *(     -1)     *aLon[i-1][j-1] +
	             ((SS*(i+1)+bottomlat)-lat)    *(1      )     *aLon[i-1][j+2] + 
	             (lat-(SS*i+bottomlat))        *(1      )     *aLon[i+2][j+2] + 
       	             (lat-(SS*i+bottomlat))        *(     -1)     *aLon[i+1][j-1] + 
	             ((SS*(i+1)+bottomlat)-lat)    *(     -1)     *aLon[i  ][j-1] +
	             ((SS*(i+2)+bottomlat)-lat)    *(     -1)     *aLon[i-1][j  ] +
	             ((SS*(i+2)+bottomlat)-lat)    *(1      )     *aLon[i-1][j+1] + 
	             ((SS*(i+1)+bottomlat)-lat)    *(1      )     *aLon[i  ][j+2] + 
	             (lat-(SS*i+bottomlat))        *(1      )     *aLon[i+1][j+2] + 
	             (lat-(SS*(i-1)+bottomlat))    *(1      )     *aLon[i+2][j+1] + 
	             (lat-(SS*(i-1)+bottomlat))    *(     -1)     *aLon[i+2][j  ];

		d[1] /= 16*SS*SS;

		if (h = al->hes) {
			h[0] =  0;
			h[1] = aLon[i+1][j+1]
			      -aLon[i  ][j+1]
			      -aLon[i+1][j  ]
			      +aLon[i  ][j  ]
			      -aLon[i+2][j-1]
			      +aLon[i-1][j-1]
			      -aLon[i-1][j+2]
			      +aLon[i+2][j+2]
			      -aLon[i+1][j-1]
			      +aLon[i  ][j-1]
			      +aLon[i-1][j  ]
			      -aLon[i-1][j+1]
			      -aLon[i  ][j+2]
			      +aLon[i+1][j+2]
			      +aLon[i+2][j+1]
			      -aLon[i+2][j  ];

			h[1] /= 16*SS*SS;
			h[2] =  0;
			}
		}
	return (rv);
	}

 void
funcadd(AmplExports *ae){
	/* Insert calls on addfunc here... */

/* Arg 3, called type, must satisfy 0 <= type <= 6:
 * type&1 == 0:	0,2,4,6	==> force all arguments to be numeric.
 * type&1 == 1:	1,3,5	==> pass both symbolic and numeric arguments.
 * type&6 == 0:	0,1	==> the function is real valued.
 * type&6 == 2:	2,3	==> the function is char * valued; static storage
			    suffices: AMPL copies the return value.
 * type&6 == 4:	4,5	==> the function is random (real valued).
 * type&6 == 6: 6	==> random, real valued, pass nargs real args,
 *				0 <= nargs <= 2.
 *
 * Arg 4, called nargs, is interpretted as follows:
 *	>=  0 ==> the function has exactly nargs arguments
 *	<= -1 ==> the function has >= -(nargs+1) arguments.
 */

	addfunc("ginv",     (rfunc)ginv,     0,  1, 0);
	addfunc("sginv",    (rfunc)sginv,    2,  1, 0);
	addfunc("hypot",    (rfunc)myhypot,  0,  2, 0);
	addfunc("ncall",    (rfunc)ncall,    0,  0, 0);
	addfunc("rncall",   (rfunc)ncall,    4,  0, 0); /*could change 4 to 6*/
	addfunc("mean0",    (rfunc)mean,     0, -1, 0);
	addfunc("mean",     (rfunc)mean,     1, -1, 0);
	addfunc("mineig",   (rfunc)mineig,   1, -1, 0);
	addfunc("kth_eig",  (rfunc)kth_eig,  1, -1, 0);
	addfunc("linear",   (rfunc)linear,   1, -1, 0);
	addfunc("C",   	    (rfunc)C,        0,  8, 0);
	addfunc("S",   	    (rfunc)S,        0,  8, 0);
	addfunc("kth_diag", (rfunc)kth_diag, 1, -1, 0);
	addfunc("kth_diag_eps",(rfunc)kth_diag_eps,1, -1, 0);
	addfunc("kth_diag2",(rfunc)kth_diag2,1, -1, 0);
	addfunc("junk",     (rfunc)junk,     1, -1, 0);
	addfunc("myerf",    (rfunc)myerf,    0,  1, 0);
	addfunc("stderf",   (rfunc)stderf,   0,  1, 0);
	addfunc("mymin",    (rfunc)mymin,    1, -1, 0);
	addfunc("nonneg",   (rfunc)nonneg,   1, -1, 0);
	addfunc("kth",      (rfunc)kth,      3, -2, 0);
	addfunc("aLatf",     (rfunc)aLatf,  0,  2, 0);
	addfunc("aLonf",     (rfunc)aLonf,  0,  2, 0);
	}

double erf (double x)
{
    return x < 0.0 ? -gammp(0.5, x*x) : gammp(0.5, x*x);
}

static double gammln(double xx)
{
    double x, tmp, ser;
    static double cof[6]={76.18009173, -86.50532033, 24.01409822,
			  -1.231739516, 0.120858003e-2, -0.536382e-5};
    int j;
    x = xx-1.0;
    tmp = x+5.5;
    tmp -= (x+0.5)*log(tmp);
    ser=1.0;
    for (j=0; j<=5;j++) {
	x += 1.0;
	ser += cof[j]/x;
    }
    return -tmp+log(2.50662827465*ser);
}

static double gammp (double a, double x)
{
    double gamser, gammcf, gln;

    if (x < 0.0 || a <= 0.0) {
	/* fprintf(Stderr,"invalid arguments in gammp\n"); exit(0); */
    }
    if (x < (a+1.0)) {
	gser(&gamser,a,x,&gln);
	return gamser;
    } else {
	gcf(&gammcf,a,x,&gln);
	return 1.0-gammcf;
    }
}

#define ITMAX 100
#define EPS 3.0e-7

static void gser(double *gamser, double a, double x, double *gln)
{
    int n;
    double sum, del, ap;
    *gln=gammln(a);
    if (x <= 0.0) {
	/* if (x < 0.0) {fprintf(Stderr,"x < 0 in gser\n"); exit(0);} */
	*gamser=0.0;
	return;
    } else {
	ap = a;
	del = sum = 1.0/a;
	for (n=1; n<= ITMAX; n++) {
	    ap += 1.0;
	    del *= x/ap;
	    sum += del;
	    if (fabs(del) < fabs(sum)*EPS) {
		*gamser=sum*exp(-x+a*log(x)-(*gln));
		return;
	    }
	}
	/*
	fprintf(Stderr,"a too large or ITMAX too small in gser\n"); exit(0);
	*/
    }
}

static void gcf( double *gammcf, double a, double x, double *gln)
{
    int n;
    double gold=0.0,g,fac=1.0, b1=1.0;
    double b0=0.0, anf, ana, an, a1, a0=1.0;
    *gln = gammln(a);
    a1 = x;
    for (n=1; n<=ITMAX; n++) {
	an= (double) n;
	ana = an - a;
	a0 = (a1+a0*ana)*fac;
	b0 = (b1+b0*ana)*fac;
	anf = an*fac;
	a1 = x*a0+anf*a1;
	b1 = x*b0+anf*b1;
	if (a1) {
	    fac = 1.0/a1;
	    g = b1*fac;
	    if (fabs((g-gold)/g) < EPS) {
		*gammcf = exp(-x+a*log(x)-(*gln))*g;
		return;
	    }
	    gold = g;
	}
    }
    /*
    fprintf(Stderr,"a too large or ITMAX too small in gcf\n"); exit(0);
    */
}

/*********************************************************************/
/***    To format the wind velocity data                           ***/
/***    All Rights Reserved                                        ***/
/*********************************************************************/

#define pi         (4*atan(1)) 
#define RADIUS     4000  
#define velconv    (3.6/1.609344)

#define gridsize    0.5
#define yy          0.75 

static void alatalon(AmplExports *ae)
{
   int   ln, lt, i, j, in, LAT, LON, index;
   char  ch[4];
   char  directory[80];
   double  lon[85][141], lat[85][141], vel[85][141], dir[85][141];
   double  aLonact[85][141], aLatact[85][141];
   double  Lonftr, Latftr;
   double  extralat, extralon;
   double  sdot, beta, sinbeta, cosbeta;
   double  weight, totalweight, templon, templat, dlat, dlon;
   FILE  *f1, *f2, *f3;

   /*---- reading the longitude file -----*/
/*   strcpy(directory,getenv("AMPLFUNC"));
   printf("AMPLFUNC = %s \n", directory);
   strcat(directory,"/flight/lon2");
   */
   strcpy(directory,"lon2");
   f1 = fopen(directory, "r");
   if (f1==NULL) printf("can not open %s\n", directory); fflush(stdout);

   for (i=0; i<85; i++) {
      fscanf(f1,"%s %d ", &ch, &in); 
      for (j=0; j<141; j++) { 
         fscanf(f1, "%lf", &lon[i][j]);
         lon[i][j] /= 10*pi/180;  /* converting into degrees */
         }
      }

   /*---- reading the latitude file ------*/
   /*
   strcpy(directory,getenv("AMPLFUNC"));
   strcat(directory,"/flight/lat2");
   */
   strcpy(directory,"lat2");
   f1 = fopen(directory, "r");
   if (f1==NULL) printf("can not open %s\n", directory); fflush(stdout);
   for (i=0; i<85; i++) {
      fscanf(f1,"%s %d", &ch, &in); 
      for (j=0; j<141; j++) {
         fscanf(f1, "%lf", &lat[i][j]);
         lat[i][j] /= 100*pi/180; /* converting into degrees */
         }
      }

   /*---- reading the magnitudes of the velocities ------*/
   /*
   strcpy(directory,getenv("AMPLFUNC"));
   strcat(directory,"/flight/vel2");
   */
   strcpy(directory,"vel2");
   f1 = fopen(directory, "r");
   if (f1==NULL) printf("can not open %s\n", directory); fflush(stdout);
   for (i=0; i<85; i++) {
      fscanf(f1,"%s %d", &ch, &in); 
      for (j=0; j<141; j++) fscanf(f1, "%lf", &vel[i][j]);    /* velocity units: m/s */
      }

   /*---- reading the directions of the velocities ------*/
   /*
   strcpy(directory,getenv("AMPLFUNC"));
   strcat(directory,"/flight/dir2");
   */
   strcpy(directory,"dir2");
   f1 = fopen(directory, "r");
   if (f1==NULL) printf("can not open %s\n", directory); fflush(stdout);
   for (i=0; i<85; i++) {
      fscanf(f1,"%s %d", &ch, &in); 
      for (j=0; j<141; j++) fscanf(f1, "%lf", &dir[i][j]); /* units: degrees */
      }

   fclose(f1);

   /*---- calculating the latitude component of the velocities ------*/
   for (i=0; i<85; i++) 
   for (j=0; j<141; j++) 
      aLatact[i][j] = sin((pi/180)*dir[i][j])*vel[i][j]*velconv/RADIUS; 

   /*---- calculating the longitude component of the velocities ------*/
   for (i=0; i<85; i++) 
   for (j=0; j<141; j++) 
      aLonact[i][j] = cos((pi/180)*dir[i][j])*vel[i][j]*velconv/(RADIUS*cos((pi/180)*lat[i][j])); 
 
   /*---- allocating memory for aLat and aLon and initializing to 0 ----*/
   LON = (int)((rightlon- leftlon)/gridsize) + 1;
   LAT = (int)((toplat- bottomlat)/gridsize) + 1;

   MALLOC(aLat, LAT, double *);
   for (lt=0; lt<LAT; lt++)   MALLOC(aLat[lt], LON, double);

   MALLOC(aLon, LAT, double *);
   for (lt=0; lt<LAT; lt++)   MALLOC(aLon[lt], LON, double);

   for (lt=0; lt< LAT; lt++) 
   for (ln=0; ln< LON; ln++) {
      aLat[lt][ln] = 0.0;   
      aLon[lt][ln] = 0.0;   
      }

   /*---- Finding the components of the velocities on the grid ------*/
   for (lt=0; lt< LAT; lt++) 
   for (ln=0; ln< LON; ln++) {   
      templon = leftlon + gridsize*ln;
      templat = bottomlat + gridsize*lt;

      totalweight = 0.0;
      for (i=0; i<85; i++)
      for (j=0; j<141; j++) 
      if (((dlat = lat[i][j] - templat) < 1) && (templat - lat[i][j] < 1)) 
      if (((dlon = lon[i][j] - templon) < 1) && (templon - lon[i][j] < 1)) {
         if (dlat < 0) dlat = -1 * dlat; 
         if (dlon < 0) dlon = -1 * dlon; 
         weight         = 1/((dlat+0.1)*(dlon+0.1)); 
         aLat[lt][ln] += aLatact[i][j] * weight;
         aLon[lt][ln] += aLonact[i][j] * weight;
         totalweight   += weight;
         }

      if (totalweight == 0.0) 
      for (i=0; i<85; i++)
      for (j=0; j<141; j++) 
      if (((dlat = lat[i][j] - templat) < 2) && (templat - lat[i][j] < 2)) 
      if (((dlon = lon[i][j] - templon) < 2) && (templon - lon[i][j] < 2)) {
         if (dlat < 0) dlat = -1 * dlat; 
         if (dlon < 0) dlon = -1 * dlon; 
         weight         = 1/((dlat+0.1)*(dlon+0.1)); 
         aLat[lt][ln] += aLatact[i][j] * weight;
         aLon[lt][ln] += aLonact[i][j] * weight;
         totalweight   += weight;
         }

      if (totalweight == 0.0) printf("%lf  %lf\n", templat, templon);
      else {
         aLat[lt][ln] /= totalweight;
         aLon[lt][ln] /= totalweight;
         }
      aLat[lt][ln] *= 180/pi;    /* converting it into degrees/hour */
      aLon[lt][ln] *= 180/pi;
      }
   
   /* to plot the velocities and directions in a vrml file */

   /* finding the future positions of the wind starting at the grid points */
   /*
   f1 = fopen("vrml.txt", "w");
   if (f1==NULL) printf("can not open vrml.txt\n"); fflush(stdout);
   f2 = fopen("index.txt", "w");
   if (f2==NULL) printf("can not open index.txt\n"); fflush(stdout);
   f3 = fopen("color.txt", "w");
   if (f3==NULL) printf("can not open color.txt\n"); fflush(stdout);

   index = 0;

   for (lt=0; lt< LAT; lt+= 2) 
   for (ln=0; ln< LON; ln+= 2) {
      templon = leftlon + gridsize*ln;    
      templat = bottomlat + gridsize*lt;

      Latftr  = 0.0;
      Lonftr  = 0.0;

      sdot = RADIUS*sqrt((pow(aLon[lt][ln]*pi/180,2)*pow(cos(templat*pi/180),2))+ pow(aLat[lt][ln]*pi/180,2));  

      if (sdot == 0) printf("1/0 in the vermal part of funcadd\n");

      sinbeta = RADIUS*aLat[lt][ln]*(pi/180)/sdot; 
      cosbeta = RADIUS*aLon[lt][ln]*(pi/180)*cos(templat*pi/180)/sdot; 

      if (cosbeta < 0)      beta = pi - asin(sinbeta);
      else if (sinbeta > 0) beta = asin(sinbeta);
      else                  beta = 2*pi + asin(sinbeta);

      Latftr = templat +(180/pi)*yy*sdot*sin(beta)/RADIUS;
      Lonftr = templon +(180/pi)*yy*sdot*cos(beta)/(RADIUS*cos(templat*pi/180));

      extralat = templat + (180/pi)*0.3*yy*sdot*sin(beta)/RADIUS;
      extralon = templon + (180/pi)*0.3*yy*sdot*cos(beta)/(RADIUS*cos(templat*pi/180));

      fprintf(f1,"%0.8f %0.8f %0.8f,\n", 2.02*sin(templon*pi/180)*cos(templat*pi/180), 2.02* sin(templat*pi/180), 2.02*cos(templon*pi/180)*cos(templat*pi/180));

      fprintf(f1,"%0.8f %0.8f %0.8f,\n", 2.02*sin(Lonftr*pi/180)*cos(Latftr*pi/180), 2.02*sin(Latftr*pi/180), 2.02*cos(Lonftr*pi/180)*cos(Latftr*pi/180)); 

      fprintf(f1,"%0.8f %0.8f %0.8f,\n", 2.02*sin(extralon*pi/180)*cos(extralat*pi/180), 2.02*sin(extralat*pi/180), 2.02*cos(extralon*pi/180)*cos(extralat*pi/180)); 

      fprintf(f2, "%d %d -1\n", index, index+1);
      fprintf(f2, "%d %d -1\n", index, index+2);

      switch ((int)(sdot*4/120)) {
         case 0: fprintf(f3, "5 0 \n"); break;
         case 1: fprintf(f3, "5 0 \n"); break;
         case 2: fprintf(f3, "5 0 \n"); break;
         case 3: fprintf(f3, "5 0 \n"); break;
         case 4: fprintf(f3, "5 0 \n"); break;
         }
      index +=3;
      }

fclose(f1);
fclose(f2);
fclose(f3);
   */
}

#ifdef __cplusplus
	}
#endif