/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:14 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_sbacc s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#include "p_sbacc.h"
/*     program DRSBACC
 *>> 1996-05-28 DRSBACC  Krogh Moved formats up.
 *>> 1994-10-19 DRSBACC  Krogh  Changes to use M77CON
 *>> 1987-12-09 DRSBACC  Lawson  Initial Code.
 *--S replaces "?": DR?BACC, ?BACC, ?BSOL
 *     Demonstration driver for SBACC & SBSOL
 *     C. L. Lawson & S. Y. Chiu, JPL, July 1987, Sept 1987.
 *     ------------------------------------------------------------------ */
		/* PARAMETER translations */
#define	C45	45.0e0
#define	HALF	0.5e0
#define	ISCALE	6
#define	LDG	24
#define	NB	2
#define	NX	10
#define	ONE	1.0e0
#define	XINC	1.0e0
#define	XLIMIT	89.5e0
#define	ZERO	0.0e0
		/* end of PARAMETER translations */
 
 
int main( )
{
	long int i, ierr2, ierr3, ig, ir, j, jt, jtprev, mt, mtotal;
	float c[NX][NX], delx, dof, dtor, g[3][LDG], rdummy, rnorm, sigfac,
	 vfac, x, yf, yfit[NX], ytrue;
	static float xtab[NX]={0.0e0,10.0e0,20.0e0,30.0e0,40.0e0,50.0e0,
	 60.0e0,70.0e0,80.0e0,90.0e0};
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const Xtab = &xtab[0] - 1;
	float *const Yfit = &yfit[0] - 1;
		/* end of OFFSET VECTORS */
 
 
	/*     ------------------------------------------------------------------ */
	printf("         Demonstration driver for SBACC and SBSOL.\n\n "
	   "Compute least-squares fit of a continuous piecewise linear\n "
	   "     function to the sine function on 0 to 90 degrees.\n \n");
	dtor = atanf( ONE )/C45;
	mtotal = 0;
	ir = 1;
	mt = 0;
	jt = 1;
	ig = 0;
	x = -XINC;
	delx = Xtab[jt + 1] - Xtab[jt];
 
L_20:
	if (x < XLIMIT)
	{
		x += XINC;
		mtotal += 1;
		if (x > Xtab[jt + 1])
		{
			printf(" Calling SBACC with JT =%3ld, MT =%3ld\n", jt, mt);
			sbacc( (float*)g, LDG, NB, &ir, mt, jt, &jtprev, &ierr2 );
			ig = ir - 1;
			mt = 0;
			jt = min( jt + 1, NX - 1 );
			delx = Xtab[jt + 1] - Xtab[jt];
		}
		ig += 1;
		mt += 1;
		g[0][ig - 1] = (Xtab[jt + 1] - x)/delx;
		g[1][ig - 1] = (x - Xtab[jt])/delx;
		g[2][ig - 1] = sinf( x*dtor );
		goto L_20;
	}
 
	if (mt > 0)
	{
		printf(" Calling SBACC with JT =%3ld, MT =%3ld\n", jt, mt);
		sbacc( (float*)g, LDG, NB, &ir, mt, jt, &jtprev, &ierr2 );
	}
 
	sbsol( 1, (float*)g, LDG, NB, ir, jtprev, yfit, NX, &rnorm, &ierr3 );
	/*                       The following statement does a type conversion. */
	dof = mtotal - NX;
	sigfac = rnorm/sqrtf( dof );
	printf(" \n MTOTAL =%4ld,  RNORM =%10.5f,  SIGFAC =%10.5f\n", mtotal, rnorm, sigfac);
	printf(" \n     X       Y        YFIT   R=Y-YFIT\n \n");
	for (i = 1; i <= NX; i++)
	{
		ytrue = sinf( Xtab[i]*dtor );
		printf(" %6.1f%10.5f%10.5f%10.5f\n", Xtab[i], ytrue, Yfit[i], Yfit[i] - ytrue);
		if (i != NX)
		{
			x = HALF*(Xtab[i + 1] + Xtab[i]);
			yf = HALF*(Yfit[i + 1] + Yfit[i]);
			ytrue = sinf( x*dtor );
			printf(" %6.1f%10.5f%10.5f%10.5f\n", x, ytrue, yf, yf - ytrue);
		}
	}
 
	/*     Compute unscaled covariance matrix in C(,).
	 * */
	for (j = 1; j <= NX; j++)
	{
		for (i = 1; i <= NX; i++)
		{
			c[j - 1][i - 1] = ZERO;
		}
		c[j - 1][j - 1] = ONE;
		sbsol( 2, (float*)g, LDG, NB, ir, jtprev, &c[j - 1][0], NX,
		 &rdummy, &ierr3 );
		sbsol( 3, (float*)g, LDG, NB, ir, jtprev, &c[j - 1][0], NX,
		 &rdummy, &ierr3 );
	}
 
	printf(" \n          Covariance matrix scaled up by 10**%1d\n \n", ISCALE);
 
	vfac = (ipow(10,ISCALE))*SQ(sigfac);
	for (i = 1; i <= NX; i++)
	{
		printf(" ");
		for (j = 1; j <= NX; j++)
		{
			printf("%7.2f", vfac*c[j - 1][i - 1]);
		}
		printf(" \n");
	}
 
	exit(0);
} /* end of function */