#include	<stdio.h>
#include	<math.h>
#define	MAXDATA	101
#define	JIKU	4
#define DATA1	7
#define	DATA	21
#define	DATA2	31

char InFile[256], OutFile[256];
FILE *in, *out;
int SilentMode = 0;

double	x[DATA], sx[DATA];
double	ai[JIKU], dai[JIKU], ao[JIKU], dao[JIKU], cao[JIKU], dcao[JIKU], agi[JIKU], dagi[JIKU], ago[JIKU], dago[JIKU];
double	cai[DATA1], dcai[DATA1], sai[DATA1], dsai[DATA1], sao[DATA1], dsao[DATA1];
double	vi, dvi, vo, dvo;
float	p[DATA], sigp[DATA];
float	bo[MAXDATA], wgt[MAXDATA], xlsq[DATA][MAXDATA], bc[MAXDATA], bd[MAXDATA], sigb[MAXDATA], bw[MAXDATA], fh[JIKU][MAXDATA];
int	iinum[MAXDATA], ih1[MAXDATA], ih2[MAXDATA], ih3[MAXDATA];

float	am[DATA2][DATA2], aam[DATA2][DATA2];

int main(int argn, char *argv[])
{
	char line[1024];

int	le[DATA1], ih[JIKU], xh[DATA1];
float	aaih[JIKU], aai[DATA1], daai[DATA1];
char	title[256];
void	sterr1(), sterr2();
double	cosaif(), sigmaf(), acosd(), trigcf(), trigaf(), trigvf();
int	i, ie, io, ip, iw, j, jj, job, jw, k, ls, n, no, np, num;
float	abc, def, ghi, pqr, xyz, bas, wl, wls, wwl, radius, sigg, zo;
float	dobs, dcal, thecc, thecl , theta;

int	m,l;

	if(argn < 3) {
		fprintf(stderr, "Usage: K.exe InputFile OutputFile [--silent]\n");
		return 1;
	}
	strcpy(InFile , argv[1]);
	strcpy(OutFile, argv[2]);
	if(argn >= 4 && strstr(argv[3], "silent")) {
		SilentMode = 1;
	}
	if((in = fopen(InFile, "r")) == NULL) {
		sprintf(line, "Can not read %s\n", InFile);
		fprintf(stderr, line);
		return 2;
	}
	if((out = fopen(OutFile, "w")) == NULL) {
		sprintf(line, "Can not write to %s\n", OutFile);
		fprintf(stderr, line);
		return 2;
	}

	if(!SilentMode) {
		fprintf(stderr, "Least Square Calculation for Lattice Constant\n");
		fprintf(stderr, "           modified by T.Kamiya    1989.10\n");
		fprintf(stderr, "           modified by T.Kamiya    2005.08.02\n\n");
	}
	fprintf(out, "Least Square Calculation for Lattice Constant\n");
	fprintf(out, "           modified by T.Kamiya    1989.10\n");
	fprintf(out, "           modified by T.Kamiya    2005.08.02\n\n");

do {
	fgets(title, 256, in);
	fprintf(out, "\n %s\n",title);

	fscanf(in, "%d %d %d %d %d %d %d %d %d",&ls,&le[1],&le[2],&le[3],&le[4],&le[5],&iw,&io,&ip);
	fprintf(out, "     HLS   LE(1)  LE(2)  LE(3)  LE(4)  LE(5)   IW     IO     IP\n");
	fprintf(out, "%7d",ls);
	for (i = 1 ; i <=5 ; i++) fprintf(out, "%7d",le[i]);
	fprintf(out, "%7d%7d%7d\n",iw,io,ip);

	num = 0;
	jj = 0;
	jw = 0;

	no = 1;
	fscanf(in, "%f%f",&wl,&radius);
	fprintf(out, "Wave Length =%f        Radius =%f\n",wl,radius);
	if (jw == 0) wwl = wl;
	jw++;
	wls = 4.0 / (wl * wl);
	fscanf(in, "%d %d %d %f",&ih[1],&ih[2],&ih[3],&bo[no]);
	while (ih[1] < 1000 || ih[1] == 2000) {
		if (ih[1] == 2000) {
			fscanf(in, "%f%f",&wl,&radius);
			fprintf(out, "Wave Length =%f        Radius =%f\n",wl,radius);
			if (jw == 0) wwl = wl;
			jw++;
			wls = 4.0 / (wl * wl);
		}
		else {
			num++;
			jj++;
			iinum[jj] = num;
			ih1[jj] =ih[1];
			ih2[jj] =ih[2];
			ih3[jj] =ih[3];
			bw[jj] = bo[no];
			switch (io) {
				case 1:	bo[no] = sin(0.01745329 * bo[no]);
					break;
				case 2:	break;
				case 3:	fscanf(in, "%f",&zo);
					xyz = cos(bo[no] / radius);
					abc = atan(zo /radius);
					def = cos(abc);
					ghi = xyz * def;
					pqr = acos(ghi);
					bo[no] = sin(0.5 * pqr);
					break;
				case 4:
					bo[no] = 0.5 * bo[no];
					bo[no] = sin(0.01745329 * bo[no]);
					break;
				case 5:
					bo[no] = wl * 0.5 / bo[no];
					break;
				default:abc=bo[no]/radius;
					xyz=1.0+abc*abc;
					bo[no]=abc/sqrt(xyz);
					break;
			}
			bas = bo[no] * bo[no];
			switch (iw) {
				case 0: fscanf(in, "%f",&sigg);sigg*=0.25;
					if (sigg <= 0) sigg = 0.25;
					wgt[no] = sigg / (bas * (1.0 - bas));
					break;
				case 1:	fscanf(in, "%f",&sigg);sigg*=2.0;
					wgt[no] = 0.25 / (bas * (1.0 - bas) * sigg * sigg);
					break;
				case 2:	wgt[no] = 1.0;
					break;
				case 3:	fscanf(in, "%f",&sigg);
					if (sigg > 0) wgt[no] = sigg;
					else wgt[no]=1.0;
					break;
				case 4:	fscanf(in, "%f",&sigg);
					wgt[no] = 1.0 / (sigg * sigg);
					break;
			}
			for (i = 1 ; i <= 3 ; i++) {
				aaih[i] = ih[i];
				fh[i][no] = aaih[i];
				xh[i] = fh[i][no] * fh[i][no];
			}
			xh[4] = 2.0 * fh[2][no] * fh[3][no];
			xh[5] = 2.0 * fh[3][no] * fh[1][no];
			xh[6] = 2.0 * fh[1][no] * fh[2][no];
			switch (ls) {
				case 1:	np = 6;
					for (i = 1 ; i <= 6 ; i++) xlsq[i][no] = xh[i];
					break;
				case 2:	np = 4;
					for (i = 1 ; i <= 3 ; i++) xlsq[i][no] = xh[i];
					xlsq[4][no] = xh[5];
					break;
				case 3:	np = 4;
					for (i = 1 ; i <= 3 ; i++) xlsq[i][no] = xh[i];
					xlsq[4][no] = xh[6];
					break;
				case 4:	np = 3;
					for (i = 1 ; i <= 3 ; i++) xlsq[i][no] = xh[i];
					break;
				case 5:	np = 2;
					xlsq[1][no] = xh[1] + xh[2];
					xlsq[2][no] = xh[3];
					break;
				case 6:	np = 1;
					xlsq[1][no] = xh[1] + xh[2] + xh[3];
					break;
				case 7:	np = 2;
					xlsq[1][no] = xh[1] + xh[2] + xh[3];
					xlsq[2][no] = xh[4] + xh[5] + xh[6];
					break;
				case 8:	np = 2;
					xlsq[1][no] = xh[1] + xh[2] + 0.5 * xh[6];
					xlsq[2][no] = xh[3];
					break;
			}
			if (le[1]) {
				np++;
				xlsq[np][no] = (1.5707963 - asin(bo[no])) * tan(1.5707963 - asin(bo[no]));
			}
			if (le[2]) {
				np++;
				xlsq[np][no] = 4.0 * bas * (1.0 - bas);
			}
			if (le[3]) {
				np++;
				xlsq[np][no] = 4.0 * bas * (1.0 - bas) * (1.0 / bo[no] + 1.0 / asin(bo[no]));
			}
			if (le[4]) {
				np++;
                                xlsq[np][no] = sin(2.0*asin(bo[no]));      
			}
			if (le[5]) {
				np++;
                                theta = asin(bo[no]);
                                xlsq[np][no] = sin(2.0*theta)*cos(theta);
			}
			bo[no] = bas * wls;
			no++;
		}
                if(feof(stdin)!=0){
			fprintf(out,"\n\nAbnormal Program END for EOF\n");
			exit(-2);
			}
		fscanf(in, "%d %d %d %f",&ih[1],&ih[2],&ih[3],&bo[no]);
	}

	--no;
	no = lsq(no, np, wl);
	fprintf(out, "             P(I)                SIGP(I)\n");
	for (i = 1 ; i <= np ; i++) fprintf(out, "          %10.7f          %10.7f\n",p[i],sigp[i]);
	fprintf(out, "         H    K    L   Dobs.   Dcal.    BO      BC       D(B)      d(2Q)\n");
	for (j = 1 ; j <= no ; j++) {
		dobs = 1.0 / (sqrt(bo[j]));
		dcal = 1.0 / (sqrt(bc[j]));
		thecc = wwl / (2.0 * dcal);
		thecl = 2.0 * asin(thecc) * 57.29578;
		fprintf(out, " (%2d)%5d%5d%5d%8.4f%8.4f%8.3f%8.3f  %9.5f  %9.5f\n",
					iinum[j],ih1[j],ih2[j],ih3[j],dobs,dcal,bw[j],thecl,bd[j],bw[j]-thecl);
	}
	for (i = 1 ; i <=3 ; i++) {
		cai[i] = 0.0;
		cao[i] = 0.0;
		dcai[i] = 0.0;
		dcao[i] = 0.0;
		agi[i] = 90.0;
		ago[i] = 90.0;
		dagi[i] = 0.0;
		dago[i] = 0.0;
	}
	switch (ls) {
		case 1:	for(i = 1 ; i <= 3 ; i++) {
				ai[i] = sqrt((double)p[i]);
				dai[i] = sigp[i] / (2.0 * ai[i]);
				aai[i] = ai[i];
				aai[i + 3] = ai[i];
				daai[i] = dai[i];
				daai[i + 3] = dai[i];
			}
			n = 3;
			for(i = 1 ; i <= 3 ; i++) {
				x[1] = p[i + 3];
				sx[1] = sigp[i + 3];
				for(j = 1 ; j <= 2 ; j++) {
					k = i + j;
					x[j + 1] = aai[k];
					sx[j + 1] = daai[k];
				}
				dcai[i] = sigmaf(&(cai[i]), n, cosaif);
			}
			recip();
			break;
		case 2:	p[5] = p[4];
			sigp[5] = sigp[4];
			p[4] = 0.0;
			sigp[4] = 0.0;
			p[6] = 0.0;
			sigp[6] = 0.0;
			for(i = 1 ; i <= 3 ; i++) {
				ai[i] = sqrt((double)p[i]);
				dai[i] = sigp[i] / (2.0 * ai[i]);
				aai[i] = ai[i];
				aai[i + 3] = ai[i];
				daai[i] = dai[i];
				daai[i + 3] = dai[i];
			}
			n = 3;
			for(i = 1 ; i <= 3 ; i++) {
				x[1] = p[i + 3];
				sx[1] = sigp[i + 3];
				for(j = 1 ; j <= 2 ; j++) {
					k = i + j;
					x[j + 1] = aai[k];
					sx[j + 1] = daai[k];
				}
				dcai[i] = sigmaf(&(cai[i]), n, cosaif);
			}
			recip();
			break;
		case 3:	p[6] = p[4];
			sigp[6] = sigp[4];
			p[4] = 0.0;
			sigp[4] = 0.0;
			p[5] = 0.0;
			sigp[5] = 0.0;
			for(i = 1 ; i <= 3 ; i++) {
				ai[i] = sqrt((double)p[i]);
				dai[i] = sigp[i] / (2.0 * ai[i]);
				aai[i] = ai[i];
				aai[i + 3] = ai[i];
				daai[i] = dai[i];
				daai[i + 3] = dai[i];
			}
			n = 3;
			for(i = 1 ; i <= 3 ; i++) {
				x[1] = p[i + 3];
				sx[1] = sigp[i + 3];
				for(j = 1 ; j <= 2 ; j++) {
					k = i + j;
					x[j + 1] = aai[k];
					sx[j + 1] = daai[k];
				}
				dcai[i] = sigmaf(&(cai[i]), n, cosaif);
			}
			recip();
			break;
		case 4:	for(i = 1 ; i <= 3 ; i++) {
				ai[i] = sqrt((double)p[i]);
				dai[i] = sigp[i] / (2.0 * ai[i]);
				ao[i] = 1.0 / ai[i];
				dao[i] = dai[i] * ao[i] * ao[i];
			}
			vi = ai[1] * ai[2] * ai[3];
			dvi = vi * sqrt((dai[1] / ai[1]) * (dai[1] / ai[1]) + (dai[2] / ai[2]) * (dai[2] / ai[2]) + (dai[3] / ai[3]) * (dai[3] / ai[3]));
			vo = 1.0 / vi;
			dvo = dvi * vo * vo;
			break;
		case 5:	ai[3] = sqrt((double)p[2]);
			dai[3] = 0.5 * sigp[2] / ai[3];
			ai[1] = sqrt((double)p[1]);
			dai[1] = 0.5 * sigp[1] / ai[1];
			ai[2] = ai[1];
			dai[2] = dai[1];
			for(j = 1 ; j <= 3 ; j++) {
				ao[j] = 1.0 / ai[j];
				dao[j] = dai[j] * ao[j] * ao[j];
			}
			vi = ai[3] * ai[1] * ai[1];
			dvi = vi * sqrt(2.0 * (dai[1] / ai[1]) * (dai[1] / ai[1]) + (dai[3] / ai[3]) * (dai[3] / ai[3]));
			vo = 1.0 / vi;
			dvo = dvi * vo * vo;
			break;
		case 6:	ai[1] = sqrt((double)p[1]);
			dai[1] = 0.5 * sigp[1] / ai[1];
			ao[1] = 1.0 / ai[1];
			dao[1] = dai[1] * ao[1] * ao[1];
			for(j = 2 ; j <= 3 ; j++) {
				ai[j] = ai[1];
				dai[j] = dai[1];
				ao[j] = ao[1];
				dao[j] = dao[1];
			}
			vi = ai[1] * ai[1] * ai[1];
			dvi = 3.0 * ai[1] * ai[1] * dai[1];
			vo = 1.0 / vi;
			dvo = dvi * vo * vo;
			break;
		case 7:	for(j = 1 ; j <= 3 ; j++) {
				ai[j] = sqrt((double)p[1]);
				cai[j] = p[2] / p[1];
				dai[j] = 0.5 * sigp[1] / ai[1];
				dcai[j] = cai[j] * sqrt((double)(sigp[1] / p[1]) * (sigp[1] / p[1]) + (sigp[2] / p[2]) * (sigp[2] / p[2]));
			}
			x[1] = ai[1];
			x[2] = cai[1];
			sx[1] = dai[1];
			sx[2] = dcai[1];
			n = 2;
			dao[1] = sigmaf(&(ao[1]), n, trigaf);
			dvi = sigmaf(&vi, n, trigvf);
			x[1] = cai[1];
			sx[1] = dcai[1];
			dcao[1] = sigmaf(&(cao[1]), n, trigcf);
			dagi[1] = sigmaf(&(agi[1]), n, acosd);
			x[1] = cao[1];
			sx[1] = dcao[1];
			dago[1] = sigmaf(&(ago[1]), n, acosd);
			for(j = 2 ; j <= 3 ; j++) {
				cao[j] = cao[1];
				ago[j] = ago[1];
				agi[j] = agi[1];
				dagi[j] = dagi[1];
				dcao[j] = dcao[1];
				dago[j] = dago[1];
				dao[j] = dao[1];
				ao[j] = ao[1];
			}
			vo = 1.0 / vi;
			dvo = dvi * vo * vo;
			break;
		case 8:	ai[1] = sqrt((double)p[1]);
			ai[2] = ai[1];
			ai[3] = sqrt((double)p[2]);
			dai[1] = 0.5 * sigp[1] / ai[1];
			dai[2] = dai[1];
			dai[3] = 0.5 * sigp[2] / ai[3];
			vi = 0.8660254 * ai[1] * ai[2] * ai[3];
			dvi = vi * sqrt(4.0 * (dai[1] / ai[1]) * (dai[1] / ai[1]) + (dai[3] / ai[3]) * (dai[3] / ai[3]));
			ao[1] = 1.1547005 / ai[1];
			ao[2] = ao[1];
			ao[3] = 1.0 / ai[3];
			dao[1] = dai[1] * ao[1] / ai[1];
			dao[2] = dao[1];
			dao[3] = dai[3] * ao[3] / ai[3];
			cai[3] = 0.5;
			cao[3] = -0.5;
			agi[3] = 60.0;
			ago[3] = 120.0;
			vo = 1.0 / vi;
			dvo = dvi * vo * vo;
			break;
	}
	fprintf(out, "Reciprocal cell constant\n");

	fprintf(out, "     a*                        b*                      c*\n");
	for(j=1;j<=3;j++) fprintf(out, "%10.7f(%10.7f)   ",ai[j],dai[j]);
	fprintf(out, "\n   alpha*                     beta*                  gamma*\n");
	for(j=1;j<=3;j++) fprintf(out, "%10.6f(%10.6f)   ",agi[j],dagi[j]);
	fprintf(out, "\n   cos a*                    cos b*                  cos c*\n");
	for(j=1;j<=3;j++) fprintf(out, "%10.6f(%10.6f)   ",cai[j],dcai[j]);
	fprintf(out, "\n   V*=%g(%g)\n",vi,dvi);
	fprintf(out, "\nDirect cell constant\n");
	fprintf(out, "       a                         b                        c\n");
	for(j=1;j<=3;j++) fprintf(out, "%8.5f(%8.5f)       ",ao[j],dao[j]);
	fprintf(out, "\n      alpha                     beta                    gamma\n");
	for(j=1;j<=3;j++) fprintf(out, "%8.4f(%8.4f)       ",ago[j],dago[j]);
	fprintf(out, "\n      cos a                    cos b                    cos c\n");
	for(j=1;j<=3;j++) fprintf(out, "%8.4f(%8.4f)       ",cao[j],dcao[j]);
	fprintf(out, "\n   V=%g(%g)\n",vo,dvo);
	if(feof(in)!=0){
		fprintf(out, "\n\nAbnormal Program END for EOF\n");
		exit(-2);
		}
	fscanf(in, "%d",&job);
} while(job != 0);

	if(!SilentMode) fprintf(stderr,"\nLattice Constant Calculation was finished!\n");
	fprintf(out,"\nLattice Constant Calculation was finished!\n");

	exit(0);
}

lsq(no, np, wl)
int 	no, np;
float	wl;
{
extern	float	p[], sigp[];
extern	float	bo[], wgt[], xlsq[][MAXDATA], bc[], bd[], sigb[], bw[];
extern	int	iinum[], ih1[], ih2[], ih3[];
extern	float	am[][DATA2], aam[][DATA2];

float	vt[DATA], oe[MAXDATA];
float	aaa, btc, dif, ot, sig, sqrsig;
int	i, j, k, nd;

	aaa = no - np;
	for(;;) {
		if (aaa < 0) {
			fprintf(out,"Number of observation is smaller than that of parameters.\n");
			fprintf(out,"The problem can not be solved\n");
			exit(-1);
		}
		else if(aaa == 0) {
			fprintf(out, "Number of observations is equal to that of parameters.\n");
			fprintf(out, "Parameters are obtained without least square calculation.\n");
			for(i = 1 ; i <= no ; i++) {
				bc[i]=bo[i];bd[i]=0;
				for(j = 1 ; j <= no ; j++) am[i][j] = xlsq[i][j];
			}
			inv30s(am, aam, np);
			for(i = 1 ; i <= np ; i++) {
				p[i] = 0.0;
				for(k = 1 ; k <= np ; k++) p[i] += aam[k][i] * bo[k];
			}
			goto ret;
		}
		else {
			for(i = 1 ; i <= np ; i++) {
				vt[i] = 0.0;
				for(k = 1 ; k <= np ; k++) {
					am[i][k] = 0.0 ;
					for(j = 1 ; j <= no ; j++) am[i][k] += xlsq[i][j] * wgt[j] * xlsq[k][j];
				}
				for(j = 1 ; j <= no ; j++) vt[i] += xlsq[i][j] * wgt[j] * bo[j];
				if (am[i][i] < 0) {
					fprintf(out, "Normal equation contains zero diagonalelement.\n");
					fprintf(out, "The problem can not be solved.\n");
					goto ret;
				}
			}
			inv30s(am, aam, np);
			for(i = 1 ; i <= np ; i++) {
				p[i] = 0.0;
				for(k = 1 ; k <= np ; k++) p[i] += aam[k][i] * vt[k];
			}
			sig = 0.0;
			for(j = 1 ; j <= no ; j++) {
				bc[j] = 0.0;
				for(i = 1 ; i <= np ; i++) bc[j] += xlsq[i][j] * p[i];
				sig += wgt[j] * (bo[j] - bc[j]) * (bo[j] - bc[j]);
			}
			sqrsig = sqrt(sig / aaa);
			ot = 0.0;
			for(j = 1 ; j <= no ; j++) {
				sigb[j] = sqrsig / sqrt(wgt[j]);
				bd[j] = bo[j] - bc[j];
				if (fabs(bd[j]) - fabs(3.0 * sigb[j]) < 0 ) oe[j] = 0;
				else {
					oe[j] = 1.0;
					ot += 1.0;
					btc = 2.0 * asin(wl * sqrt(bc[j]) / 2.0);
					btc /= 0.0174533;
					dif = bw[j] - btc;
					fprintf(out, "This observation is eliminated.\n");
                    fprintf(out, "No. = %2d  bo = %10.3f  bc = %10.3f dif(B) = %10.3f\n", j, bw[j], btc, dif);
				}
			}
			if (ot != 0) {
				k = 1;
				for(j = 1 ; j <= no ; j++) {
					if (oe[j] == 0) {
						bo[k] = bo[j];
						wgt[k] = wgt[j];
						iinum[k] = iinum[j];
						ih1[k] = ih1[j];
						ih2[k] = ih2[j];
						ih3[k] = ih3[j];
						bw[k] = bw[j];
						for(i = 1 ; i <= np ; i++) xlsq[i][k] = xlsq[i][j];
						k++;
					}
				}
				no = k - 1;
			}
			else {
				for(i = 1 ; i <=np ; i++) sigp[i] = sqrt(aam[i][i]) * sqrsig;
				goto ret;
			}
		}
	}
	ret:	return(no);
}

double	volume()
{
extern	double	x[];

	return(x[1] * x[2] * x[3] * sqrt(1.0 - x[4] * x[4] - x[5] * x[5] - x[6] * x[6] + 2.0 * x[4] * x[5] * x[6]));
}

double	sinal()
{
extern	double	x[];

	return(sqrt(1.0 - x[1] * x[1]));
}

double	sigmaf(double *y, int n, double (*f)())
{
extern	double	x[], sx[];
double	ff,xt, sigy;
int	i;

	sigy = 0.0;
	*y = (*f)();
	for (i = 1 ; i <= n ; i++) {
		xt = x[i];
		x[i] = 0.0001 * sx[i] + x[i];
		ff=((*f)()-*y)*10000.0;
		sigy += ff*ff;
		x[i] = xt;
	}
	return(sqrt(sigy));
}

inv30s(a, b, n)
float	a[][DATA2], b[][DATA2];
int	n;
{
double	dn;
float	r,rr,s,t;
int	i,ii,j,k,l,m;

	dn = 0.000000001;
	for(m = 1 ; m <= n ; m++) {
		for (l = 1 ; l <= n ; l++) b[m][l] = 0;
	}
	for(m = 1 ; m <= n ; m++) b[m][m] = 1;
	for(i = 1 ; i <= n ; i++) {
		if ((fabs(a[i][i]) - dn ) < 0) {
			ii = i + 1;
			j = ii;
			while ((fabs(a[i][j]) - dn) < 0) {
				if (j > n) {
					fprintf (out, "No Solution !!");
					exit(-1);
				}
				else j++;
			}
			for (k = 1 ; k <= n ; k++) {
				r = a[k][i];
				a[k][i] = a[k][j];
				a[k][j] = r;
				rr = b[k][i];
				b[k][i] = b[k][j];
				b[k][j] = rr;
			}
		}
		s = a[i][i];
		for(j = 1 ; j <= n ; j++) {
			a[i][j] = a[i][j] / s;
			b[i][j] = b[i][j] / s;
		}
		for(k = 1 ; k <= n ; k++) {
			if (k != i) {
				t = a[k][i];
				for (l = 1 ; l <= n ; l++) {
					a[k][l] = a[k][l] - t * a[i][l];
					b[k][l] = b[k][l] - t * b[i][l];
				}
			}
		}
	}
	return;
}

double	acosd()
{
extern	double	x[];
double	ac;

	if (x[1] == 0) ac = 90.0;
	else ac = acos(x[1]) * 57.29578;
	return (ac);
}

double	cosaif()
{
extern	double	x[];

	return (x[1] / (x[2] * x[3]));
}

double	trigcf()
{
extern	double	x[];

	return (x[1] * (x[1] - 1.0) / (1.0 - x[1] * x[1]));
}

double	trigaf()
{
extern	double	x[];

	return (sqrt((1 - x[2] * x[2]) / (1.0 - 3.0 * x[2] * x[2] + 2.0 * x[2] * x[2] * x[2])) / x[1]);
}

double	trigvf()
{
extern	double	x[];
	return (x[1]*x[1]*x[1]*sqrt(1.0-3.0*x[2]*x[2]+2.0*x[2]*x[2]*x[2]));
}

recip()
{
extern	double	x[], sx[];
extern	double	ai[], dai[], ao[], dao[], cao[], dcao[], agi[], dagi[], ago[], dago[];
extern	double	cai[], dcai[], sai[], dsai[], sao[], dsao[];
extern	double	vi, dvi, vo, dvo;

double	rax(), rco(), volume(), sinal(), acosd();
int	j, jk, k, m, n;

	n = 1;
	for(j = 1 ; j <=3 ; j++) {
		x[1] = cai[j];
		sx[1] = dcai[j];
		dsai[j] = sigmaf(&(sai[j]), n, sinal);
		cai[j + 3] = cai[j];
		dcai[j + 3] = dcai[j];
		sai[j + 3] = sai[j];
		dsai[j + 3] = dsai[j];
	}
	m = 3;
	for(j = 1 ; j <=3 ; j++) {
		for(k = 1 ; k <=3 ; k++) {
			jk = j + k;
			x[k] = cai[jk - 1];
			sx[k] = dcai[jk - 1];
		}
		dcao[j] = sigmaf(&(cao[j]), m, rco);
		x[1] = cao[j];
		sx[1] = dcao[j];
		dsao[j] = sigmaf(&(sao[j]), n, sinal);
		sao[j + 3] = sao[j];
		dsao[j + 3] = dsao[j];
	}
	n = 4;
	for(j = 1 ; j <=3 ; j++) {
		x[1] = ai[j];
		x[2] = cai[j];
		x[3] = cai[j + 1];
		x[4] = cai[j + 2];
		sx[1] = dai[j];
		sx[2] = dcai[j];
		sx[3] = dcai[j + 1];
		sx[4] = dcai[j + 2];
		dao[j] = sigmaf(&(ao[j]), n, rax);
	}
	for(j = 1 ; j <=3 ; j++) {
		x[j] = ai[j];
		x[j + 3] = cai[j];
		sx[j] = dai[j];
		sx[j + 3] = dcai[j];
	}
	n = 6;
	dvi = sigmaf(&vi, n, volume);
	vo = 1.0 / vi;
	dvo = dvi * vo * vo;
	n = 1;
	for(j = 1 ; j <=3 ; j++) {
		x[1] = cai[j];
		sx[1] = dcai[j];
		dagi[j] = sigmaf(&(agi[j]), n, acosd);
		x[1] = cao[j];
		sx[1] = dcao[j];
		dago[j] = sigmaf(&(ago[j]), n, acosd);
	}
}

double	rax()
{
extern	double	x[];

	return(sqrt((1.0 - x[2] * x[2]) / (1.0 - x[2] * x[2] - x[3] * x[3] -x[4]* x[4] + 2.0 * x[2] * x[3] * x[4]))/x[1]);
}

double	rco()
{
extern	double	x[];

	return((x[2] * x[3] -x[1]) / sqrt((1.0 - x[2] * x[2]) * (1.0 - x[3] * x[3])));
}