/*
	xy.c

	this is	'latlon.c'  with XY defined

	x1.bat = cl /AL xy.c proj1 getfile

*/
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "proj1.h"
#include "nomouse.h"

#define XY 1

#define R      6367.65
#define PI		3.1415927
#define PI2		1.5707963
#define PI4		0.7853981
#define DEG_RAD 0.01745329251				/* converts degrees to radians */
#define KM_MILE 0.621371192				/* converts km to miles */
#define BLACK     0
#define D_GREY    1
#define L_GREY    2
#define WHITE     3
#define RED       4
#define ORANGE    5
#define YELLOW    6
#define GREEN     7
#define CYAN      8
#define BLUE      9
#define MAGENTA  10
#define TAN      11
#define BROWN    12
#define L_BLUE   13
#define BUFF     14
#define D_BLUE   15
#define MAX_PTS  30

int Imvis=0;
double PixelsPerInch;
int Times=20;
int ScreenXs=640,ScreenYs=480;
double Lat[MAX_PTS],Lon[MAX_PTS],X[MAX_PTS],Y[MAX_PTS];
double Tx[MAX_PTS],Ty[MAX_PTS];
struct control Constant;
unsigned char Buffer[12][640];
FILE *FpTxt;
char String[10][80];

extern struct projection *Proj;

int translate(struct control *cc,double *lat,double *lon,
		double *x,double *y,double *tx,double *ty,int max,FILE *fpin,
		char *file,FILE *fptxt);
double rotate(int num_pts,double *lat,double *lon,double *x,double *y,
	double *tx,double *ty,struct control *cc,double *xbar,double *ybar,
	double dscale,double *theta,int times,double xcor);
double stretch(int num_pts,double *lat,double *lon,double *x,double *y,
	double *tx,double *ty,struct control *cc,double *xbar,double *ybar,
	double dscale,double theta,int times,double *xcor);
double set_scale(int num_pts,double *lat,double *lon,double *x,double *y,
	double *tx,double *ty,struct control *cc,double *xbar,double *ybar,
	double dscale,double theta,int times,double xcor);
double get_dist(int num_pts,double *lat,double *lon,double *x,double *y,
	double *tx,double *ty,struct control *cc,double *xbar,double *ybar,
	double theta,double xcor);
int set_trig(double theta,double *s,double *sint,double *cost);
double digit_x(double mx,double my,double xb,
		double sint,double cost,double s,int rot,double xcor);
double digit_y(double mx,double my,double yb,
		double sint,double cost,double s,int rot);
double map_x(double x,double y,double xb,double yb,
		double sint,double cost,double xcor);
double map_y(double x,double y,double xb,double yb,
		double sint,double cost,double xcor);
long lvl(long attrib,long color);

int P;

main(int argc,char *argv[])
{
	int i,j,k,n;
	char ans,filename[100];
	int pnum,ret;
	double lat,lon,x,y;
	int xc=0,yc=0,xs=640,ys=480,speed=8;
	double clat,clon,cx,cy;
	FILE *fp,*fpin,*fptxt;
	char file[100],string[100];

	if(argc>1)
		strcpy(file,argv[1]);
	else
	{
		printf("Choose an '.xyc' file to be used to\n");
		printf("navigate the lat.lon points in the\n");
		printf(".'raw' file you will choose later.\n\n\n");
		if(get_file_name("*.xyc",file)<0)
		{
			printf("No files chosen.\n\n");
			exit(0);
		}
/*		printf("Give input file name.  Type will be '.xyc'.\n");
		scanf("%s",file);*/
	}
	for(k=0;k<strlen(file);k++)
		if(file[k]=='.')
			file[k]='\0';
	sprintf(string,"%s.txt",file);
	fptxt=fopen(string,"wt");
	strcpy(string,file);
	strcat(string,".xyc");
	fpin=fopen(string,"rt");
	if(!fpin)
	{
		printf("Could not open file '%s'.\n",string);
		exit(0);
	}
	pnum=choose_projection(&Constant,fptxt);
	if(pnum>=0)
	{
		printf("Constant.proj = %c\n",Constant.proj);
		printf("Give parameters.\n\n");
		get_constants(&Constant,fptxt);
		translate(&Constant,Lat,Lon,X,Y,Tx,Ty,MAX_PTS,fpin,file,fptxt);
		fclose(fptxt);
	}
}


/********************************************************************
**
**
**
********************************************************************* */

int translate(struct control *cc,double *lat,double *lon,
		double *x,double *y,double *tx,double *ty,int max,FILE *fpin,
		char *file,FILE *fptxt)

{
	int i,j,k,n;
	FILE *fpout,*fpraw,*fphdr;
	char string[100],cr=0xd,name[100],cval;
	int num_pts=-1,val,ret,rot,pval;
	long lval,color,attrib;
	double dist,xb,yb,dnum,dscale=2.0,scale,d,d_up,d_down,error,dxcor=2.0;
	double theta=0.0,dtheta=10.0*DEG_RAD,xcor=1.0;
	double llat,llon,lx,ly,mx,my,cost,sint,s;
	double xdist,ydist,dsc=R*PI/180.0;
	int meter=0,foot=0;

#ifdef XY
	fphdr=fopen("xy.hdr","wt");
	if(!fphdr)
	{
		printf("Could not write on disk.\n\n");
		exit(0);
	}
#endif
	i=0;
	while((val=fscanf(fpin,"%s",string))==1&&string[0]!='*')
	{
		if(i>=max)
		{
			printf("There are more than %d control points.\n\n",max);
			return(-3);
		}
		sscanf(string,"%lf",lat+i);
		fscanf(fpin,"%lf%lf%lf",lon+i,x+i,y+i);
		printf("%9.5lf %10.5lf %8.2lf %8.2lf\n",lat[i],lon[i],x[i],y[i]);
#ifdef XY
		fprintf(fphdr,"%9.5lf %10.5lf\n %10.5lf %10.5lf\n",
			lat[i],lon[i],x[i],y[i]);
		i+=1;
#endif
	}
#ifdef XY
	fprintf(fphdr,"%s\n",string);
	fclose(fpin);
	fclose(fphdr);
	fphdr=fopen("xy.hdr","rt");
	if(!fphdr)
	{
		printf("Error opening header file to read.\n\n");
		exit(0);
	}
#endif
	if(strlen(string)>1)
	{
		if(sscanf(string+1,"%lf",&PixelsPerInch)>0)
		{
			if(PixelsPerInch>0.0)
				Imvis=1;	/* This file was created by program 'imvis.c' */	
		}
	}
	num_pts=i;

/*	test_grat(num_pts,lat,lon,x,y,*cc);*/

	dnum=num_pts;
	if(val==0)
		return(-2);
	printf("There are %d control points.\n\n",num_pts);
	theta=0.0;
	cc->scale=1.0;
	error=set_scale(num_pts,lat,lon,x,y,tx,ty,cc,&xb,&yb,dscale,theta,1,xcor);
	printf("Error after 'set_scale' = %lf\n",error);
	for(i=0;i<Times;i++)
	{
error=set_scale(num_pts,lat,lon,x,y,tx,ty,cc,&xb,&yb,dscale,theta,1,xcor);
		printf("%cscale = %10.1lf  ",cr,cc->scale);
error=rotate(num_pts,lat,lon,x,y,tx,ty,cc,&xb,&yb,dtheta,&theta,1,xcor);
		printf("theta = %10.6lf  ",theta/DEG_RAD);
error=stretch(num_pts,lat,lon,x,y,tx,ty,cc,&xb,&yb,dxcor,theta,1,&xcor);
		printf("x stretch = %10.6lf",xcor);
		dscale=(dscale-1.0)/2.0+1.0;
		dtheta/=2.0;
		dxcor=(dxcor-1.0)/2.0+1.0;
	}
	printf("\n\nx correction = %8.6lf\n",xcor);
	printf("rotation = %12.5lf deg.\n",theta/DEG_RAD);
	if(Imvis==0)
	{
		printf("\nscale = %8.5lf km/count   or %10.4lf km/mm   or %0.0lf:1\n",
		R/cc->scale,R/cc->scale*40.0,R/cc->scale*40000000.0);
		printf("tablet root mean error = %10.3lf mm\n",error/40.0);
	}
	else
	{
		printf("\nscale = %8.5lf km/pixel   or %10.4lf km/mm   or %0.0lf:1\n",
			R/cc->scale,R/cc->scale*PixelsPerInch/25.4,
			R/cc->scale*PixelsPerInch*39370.07874);
		printf("image scale = %0.3lf Pixels/Inch\n",PixelsPerInch);
		printf("image root mean error = %10.3lf mm\n",error/PixelsPerInch*25.4);
		printf("image root mean error = %10.3lf pixels\n",error);
	}
	printf("earth root mean error = ");
	if(error*R/cc->scale>1.0)
		printf("%10.3lf km = ",error*R/cc->scale);
	else
	{
		printf("%10.3lf m = ",error*R/cc->scale*1000.0);
		meter=1;
	}
	if(error*R/cc->scale*KM_MILE>1.0)
		printf("%10.3lf mi\n",error*R/cc->scale*KM_MILE);
	else
	{
		printf("%10.3lf ft\n",error*R/cc->scale*KM_MILE*5280.0);
		foot=1;
	}
	printf("HIT ANY KEY TO GET lat/lon comparisons\n\n");getch();
	i=0;
	cost=cos(theta*(-1.0));
	sint=sin(theta*(-1.0));
	printf("         input*                output**      ");
	if(meter==0)
		printf("           errors (km)\n");
	else
		printf("           errors (m) \n");
	printf("-- ------------------   -------------------  ");
	printf(" --------------------------\n");
	printf("      Lat      Lon         Lat       Lon     ");
	printf("    X        Y       Total\n");
	printf("-- -------- ---------   --------- ---------  ");
	printf(" -------  -------   -------\n");

	if(fptxt)
		fprintf(fptxt,
			"\nx correction = %8.6lf       (xcor)\n",xcor);
	if(fptxt)
		fprintf(fptxt,
			"    rotation = %8.6lf deg.  (theta)\n",theta/DEG_RAD);
	if(Imvis==0)
	{
		if(fptxt)
			fprintf(fptxt,
			"\nscale = %8.5lf km/count   or %10.4lf km/mm   or %0.0lf:1\n",
		R/cc->scale,R/cc->scale*40.0,R/cc->scale*40000000.0);
		if(fptxt)
			fprintf(fptxt,
			"tablet root mean error = %10.3lf mm\n",error/40.0);
	}
	else
	{
		if(fptxt)
			fprintf(fptxt,
			"\nscale = %8.5lf km/pixel   or %10.4lf km/mm   or %0.0lf:1\n",
			R/cc->scale,R/cc->scale*PixelsPerInch/25.4,
			R/cc->scale*PixelsPerInch*39370.07874);
		if(fptxt)
			fprintf(fptxt,
			"image scale = %0.3lf Pixels/Inch\n",PixelsPerInch);
		if(fptxt)
			fprintf(fptxt,
			"image root mean error = %10.3lf mm\n",error/PixelsPerInch*25.4);
		if(fptxt)
			fprintf(fptxt,
			"image root mean error = %10.3lf pixels\n",error);
	}
	if(fptxt)
	{
		fprintf(fptxt,"earth root mean error = ");
		if(meter==0)
			fprintf(fptxt,"%10.3lf km = ",error*R/cc->scale);
		else
			fprintf(fptxt,"%10.3lf m = ",error*R/cc->scale*1000.0);
		if(foot==0)
			fprintf(fptxt,"%10.3lf mi\n",error*R/cc->scale*KM_MILE);
		else
			fprintf(fptxt,"%10.3lf ft\n",error*R/cc->scale*KM_MILE*5280.0);
		fprintf(fptxt,"X0 = %10.3lf  Y0 = %10.3lf\n",xb,yb);
	}
	if(fptxt)
	{
		fprintf(fptxt,"\n---------------------------------------------");
		fprintf(fptxt,"---------------------------\n");
		fprintf(fptxt,"         input*                output**      ");
		if(meter==0)
			fprintf(fptxt,"           errors (km)\n");
		else
			fprintf(fptxt,"           errors (m) \n");
		fprintf(fptxt,"-- ------------------   -------------------  ");
		fprintf(fptxt," --------------------------\n");
		fprintf(fptxt,"      Lat      Lon         Lat       Lon     ");
		fprintf(fptxt,"    X        Y       Total\n");
		fprintf(fptxt,"-- -------- ---------   --------- ---------  ");
		fprintf(fptxt," -------  -------   -------\n");
	}
	set_trig(theta,&s,&sint,&cost);
	for(i=0;i<num_pts;i++)
	{
		lx=x[i];
		ly=y[i];
		llat=llon=0.0;
		mx=map_x(x[i],y[i],xb,yb,sint,cost,xcor);
		my=map_y(x[i],y[i],xb,yb,sint,cost,xcor);
		ret=xy2ll(mx,my,&llat,&llon,*cc);
		if(ret>=0)
		{
			xdist=(llon-lon[i])*cos(lat[i]*DEG_RAD)*dsc;
			ydist=(llat-lat[i])*dsc;
			dist=sqrt(xdist*xdist+ydist*ydist);
			if(meter==0)
				fprintf(fptxt,
				"%2d %8.4lf %9.4lf -> %8.4lf %9.4lf: %8.3lf %8.3lf  %8.3lf\n",
					i+1,lat[i],lon[i],llat,llon,xdist,ydist,dist);
			else
				fprintf(fptxt,
				"%2d %8.4lf %9.4lf -> %8.4lf %9.4lf: %8.2lf %8.2lf  %8.2lf\n",
					i+1,lat[i],lon[i],llat,llon,
					xdist*1000.0,ydist*1000.0,dist*1000.0);
			if(meter==0)
				printf(
				"%2d %8.4lf %9.4lf -> %8.4lf %9.4lf: %8.3lf %8.3lf  %8.3lf\n",
					i+1,lat[i],lon[i],llat,llon,xdist,ydist,dist);
			else
				printf(
				"%2d %8.4lf %9.4lf -> %8.4lf %9.4lf: %8.2lf %8.2lf  %8.2lf\n",
					i+1,lat[i],lon[i],llat,llon,
					xdist*1000.0,ydist*1000.0,dist*1000.0);
		}
	}
	if(fptxt)
	{
		fprintf(fptxt,
		"\n---------------------------------------------");
		fprintf(fptxt,
		"---------------------------\n");
		fprintf(fptxt," *  control point latitude and longitude coordinates\n");
		fprintf(fptxt,
			"**  computed from the corresponding x and y coordinates\n\n");
		fprintf(fptxt,
			"X and Y values from the '.xyc' file are corrected before\n");
		fprintf(fptxt,
			"converting to lat/lon for rotation (theta), directional\n");
		fprintf(fptxt,
			"stretching (xcor), translation (X0 and Y0), scale and\n");
		fprintf(fptxt,
			"for the fact that Y is down int the '.xyc' file.\n\n");
		fprintf(fptxt,"   x = (X+X0)*xcor*cos(-theta)-(Y+Y0)*sin(-theta)\n");
		fprintf(fptxt,"   y =-(X+X0)*xcor*sin(-theta)-(Y+Y0)*cos(-theta)\n\n");
		fprintf(fptxt,"   X = x*cos(theta)/xcor-y*sin(theta)/xcor-X0\n");
		fprintf(fptxt,"   Y =-x*sin(theta)-y*cos(theta)-Y0\n\n");
		fprintf(fptxt,"   where:\n");
		fprintf(fptxt,"      X = column of image or distance right on tablet\n");
		fprintf(fptxt,"      Y = row of image or distance left on tablet\n");
		fprintf(fptxt,"     X0 = right offset to projected map origin\n");
		fprintf(fptxt,"     Y0 = down offset to projected map origin\n");
		fprintf(fptxt,"      x = meters left on projected map\n");
		fprintf(fptxt,"      y = meters up on projected map\n\n");
		fprintf(fptxt,"   functions used to translate systems:\n\n");
		fprintf(fptxt,
"int xy2ll(double x,double y,double *lat,double *lon,struct control cc)\n");
		fprintf(fptxt,
"int ll2xy(double lat,double lon,double *x,double *y,struct control cc)\n");
	}
	sprintf(string,"%s.raw",file);
#ifdef XY
	do
	{
		n=0;
		printf("\nChoose a '.raw' file to process.\n\n");
		if(get_file_name("*.raw",name)<0)
		{
			printf("No files chosen.\n\n");
			exit(0);
		}
/*		printf("Give input file name (type will be made .raw)\n\n");
		scanf("%s",name);*/
		for(k=0;k<strlen(name);k++)
			if(name[k]=='.')
				name[k]='\0';
		sprintf(string,"%s.raw",name);
		fpraw=fopen(string,"rt");
		if(fpraw)
		{
			printf(
"Give color (0 - 255 from image color palette) for plotting lines on image\n\n");
			printf(
"     Attributes will be automatically appended in XYC file.\n\n");
			scanf("%ld",&color);
			rot=set_trig(theta,&s,&sint,&cost);
			sprintf(string,"%s.xyc",name);
			fpout=fopen(string,"wt");
			if(!fpout)
			{
				printf("Could not open '%s' to write.\n\n",string);
				exit(0);
			}
			rewind(fphdr);
			while((cval=fgetc(fphdr))!='*'&&cval!=EOF)
				if(cval!=0xd)
					fprintf(fpout,"%c",cval);
			fprintf(fpout,"*\n");
			i=0;
			printf("--------------------------------\n");
			printf("Number Attribute Color XYC Entry\n");
			printf("------ --------- ----- ---------\n");
			while(fscanf(fpraw,"%ld",&attrib)==1)
			{
				fprintf(fpout,"%4ld",lvl(attrib,color));
				printf("%c%6d %9ld%6ld%10ld",cr,++n,attrib,color,lvl(attrib,color));
				do
				{
					if(i%3==0)
						fprintf(fpout,"\n");
					else
						fprintf(fpout,"  ");
					fscanf(fpraw,"%lf%s",&llat,string);
					sscanf(string,"%lf",&llon);
					i+=1;
					ll2xy(llat,llon,&mx,&my,*cc);
					lx=digit_x(mx,my,xb,sint,cost,s,rot,xcor);
					ly=digit_y(mx,my,yb,sint,cost,s,rot);
					fprintf(fpout,"%8.4lf %9.4lf",lx,ly);
					if(string[strlen(string)-1]=='/')
						i=0;
				}while(i>0);
				fprintf(fpout,"//\n");
			}
			printf("\n");
			fclose(fpraw);
			fclose(fpout);
		}
		else
			printf("Could not open '%s' to read.\n\n");
		printf("Do you want to process another '.raw' file?  (y or n)\n\n");
	}while(getch()!='n');
#else
	printf("\nDo you want to create '%s'?  (y or n)\n",string);
	if(getch()!='y')
	{
		fclose(fpin);
		return(-1);
	}
	fpout=fopen(string,"wt");
	if(!fpout)
	{
		printf("Could not open '%s' to write.\n\n",string);
		exit(0);
	}
	i=0;
	while(fscanf(fpin,"%d",&val)==1)
	{
		fprintf(fpout,"%4d",val);
		do
		{
			if(i%3==0)
				fprintf(fpout,"\n");
			else
				fprintf(fpout,"  ");
			fscanf(fpin,"%lf%s",&lx,string);
			sscanf(string,"%lf",&ly);
			i+=1;
			set_trig(theta,&s,&sint,&cost);
			mx=map_x(lx,ly,xb,yb,sint,cost,xcor);
			my=map_y(lx,ly,xb,yb,sint,cost,xcor);
			xy2ll(mx,my,&llat,&llon,*cc);
			fprintf(fpout,"%8.4lf %9.4lf",llat,llon);
			if(string[strlen(string)-1]=='/')
				i=0;
		}while(i>0);
		fprintf(fpout,"//\n");
	}
	fclose(fpin);

#endif	
	fclose(fpout);
}

/**********************************************************************
**
**
**
**********************************************************************/

long lvl(long attrib,long color)

{
	long lval,sign=1;

	if(attrib<0)
		sign=-1;
	lval=(attrib*sign*1000+color)*sign;
	return(lval);
}



/********************************************************************
**
**
**
********************************************************************* */

double rotate(int num_pts,double *lat,double *lon,double *x,double *y,
	double *tx,double *ty,struct control *cc,double *xbar,double *ybar,
	double dscale,double *theta,int times,double xcor)

{
	int i,j,k;
	double dist,xb,yb,dnum,dtheta=dscale*DEG_RAD,d,d_up,d_down,theta0;
	int time=0;

	for(k=0;k<times;k++)
	{
		time=0;
		do
		{
			time+=1;
			d=get_dist(num_pts,lat,lon,x,y,tx,ty,cc,&xb,&yb,*theta,xcor);
			theta0=*theta;
			*theta=theta0+dtheta;
			d_up=get_dist(num_pts,lat,lon,x,y,tx,ty,cc,&xb,&yb,*theta,xcor);
			*theta=theta0-dtheta;
			d_down=get_dist(num_pts,lat,lon,x,y,tx,ty,cc,&xb,&yb,*theta,xcor);
			*theta=theta0;
			if(d_up<d)
				*theta=theta0+dtheta;
			else if(d_down<d)
				*theta=theta0-dtheta;
		}while((d_up<d||d_down<d)&&time<100);
		dtheta/=10.0;
	}
	*xbar=xb;
	*ybar=yb;
	return(d);
}

/********************************************************************
**
**
**
********************************************************************* */

double stretch(int num_pts,double *lat,double *lon,double *x,double *y,
	double *tx,double *ty,struct control *cc,double *xbar,double *ybar,
	double dxcor,double theta,int times,double *xcor)

{
	int i,j,k;
	double dist,xb,yb,dnum,d,d_up,d_down,xcor0;
	int time=0;

	for(k=0;k<times;k++)
	{
		time=0;
		do
		{
			time+=1;;
			d=get_dist(num_pts,lat,lon,x,y,tx,ty,cc,&xb,&yb,theta,*xcor);
			xcor0=*xcor;
			*xcor=xcor0*dxcor;
			d_up=get_dist(num_pts,lat,lon,x,y,tx,ty,cc,&xb,&yb,theta,*xcor);
			*xcor=xcor0/dxcor;
			d_down=get_dist(num_pts,lat,lon,x,y,tx,ty,cc,&xb,&yb,theta,*xcor);
			*xcor=xcor0;
			if(d_up<d)
				*xcor=xcor0*dxcor;
			else if(d_down<d)
				*xcor=xcor0/dxcor;
		}while((d_up<d||d_down<d)&&time<100);
		dxcor/=10.0;
	}
	*xbar=xb;
	*ybar=yb;
	return(d);
}


/********************************************************************
**
**
**
********************************************************************* */

double set_scale(int num_pts,double *lat,double *lon,double *x,double *y,
	double *tx,double *ty,struct control *cc,double *xbar,double *ybar,
	double dscale,double theta,int times,double xcor)

{
	int i,j,k;
	double dist,xb,yb,dnum,scale,d,d_up,d_down;
	int time=0;

	for(k=0;k<times;k++)
	{
		time=0;
		do
		{
			time+=1;
			d=get_dist(num_pts,lat,lon,x,y,tx,ty,cc,&xb,&yb,theta,xcor);
			scale=cc->scale;
			cc->scale=scale*dscale;
			d_up=get_dist(num_pts,lat,lon,x,y,tx,ty,cc,&xb,&yb,theta,xcor);
			cc->scale=scale/dscale;
			d_down=get_dist(num_pts,lat,lon,x,y,tx,ty,cc,&xb,&yb,theta,xcor);
			cc->scale=scale;
			if(d_up<d)
				cc->scale=scale*dscale;
			else if(d_down<d)
				cc->scale=scale/dscale;
		}while((d_up<d||d_down<d)&&time<100);
		dscale=(dscale-1.0)/2.0+1.0;
	}
	*xbar=xb;
	*ybar=yb;
	return(d);
}

/********************************************************************
**
**
**
********************************************************************* */

double get_dist(int num_pts,double *lat,double *lon,double *x,double *y,
	double *tx,double *ty,struct control *cc,double *xbar,double *ybar,
	double theta,double xcor)

{
	int i,j,k;
	double xb=*xbar,yb=*ybar,dist,dnum=num_pts;
	double s,sint,cost,mx[MAX_PTS],my[MAX_PTS],xd,yd;
	int rot;

	xb=yb=0.0;
	rot=set_trig(theta,&s,&sint,&cost);
	for(i=0;i<num_pts;i++)
	{
		ll2xy(lat[i],lon[i],mx+i,my+i,*cc);
		tx[i]=digit_x(mx[i],my[i],(double)0.0,sint,cost,s,rot,xcor);
		ty[i]=digit_y(mx[i],my[i],(double)0.0,sint,cost,s,rot);
		xb+=(tx[i]-x[i])/dnum;
		yb+=(ty[i]-y[i])/dnum;
	}
	dist=0.0;
	for(i=0;i<num_pts;i++)
	{
		xd=digit_x(mx[i],my[i],xb,sint,cost,s,rot,xcor);
		yd=digit_y(mx[i],my[i],yb,sint,cost,s,rot);
		dist+=(x[i]-xd)*(x[i]-xd)/dnum+
				(y[i]-yd)*(y[i]-yd)/dnum;
	}
	*xbar=xb;
	*ybar=yb;
	return(sqrt(dist));
}


/**********************************************************************
**
**
**
**********************************************************************/

double map_y(double x,double y,double xb,double yb,
		double sint,double cost,double xcor)

{
	double my;

	my=-(x+xb)*xcor*sint-(y+yb)*cost;
	return(my);
}

/**********************************************************************
**
**
**
**********************************************************************/

double map_x(double x,double y,double xb,double yb,
		double sint,double cost,double xcor)

{
	double mx;

	mx=(x+xb)*xcor*cost-(y+yb)*sint;
	return(mx);
}



/**********************************************************************
**
**
**
**********************************************************************/

double digit_y(double mx,double my,double yb,
		double sint,double cost,double s,int rot)

{
	double y;

	y=-mx*sint-my*cost-yb;	
	return(y);
}

/**********************************************************************
**
**
**
**********************************************************************/

double digit_x(double mx,double my,double xb,
		double sint,double cost,double s,int rot,double xcor)

{
	double x;

	x=mx*cost/xcor-my*sint/xcor-xb;
	return(x);
}


/**********************************************************************
**
**
**
**********************************************************************/

int set_trig(double theta,double *s,double *sint,double *cost)

{
	double ct,st;

	if(theta==0.0)
	{
		*s=0.0;
		*sint=0.0;
		*cost=1.0;
		return(0);
	}
	ct=cos(theta);
	st=sin(theta);
	*s=(ct/st+st/ct);
	*sint=st;
	*cost=ct;

	return(1);
}