/*
	pjtest.c

	pjt1.bat = cl /AL pjtest.c project vimage2 graphlib

*/
#include <stdio.h>
#include <math.h>

#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 Mouse=0;
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
{
	char proj;
	int pnum;
	double lon_0,lon_1,lon_2,lon_p;
	double lat_0,lat_1,lat_2,lat_p;
	double pixsize,scale;
	double x0,y0;
	double c,C,F,H,k,k0,k1,k2,n,P,rho_0;
	double omega; /* elevation (tilt of surface plane) */
	double gamma; /* azimuth (east of north) */
	double horiz; /* radius of horizon arc */
}Constant;

struct projection
{
	char name[80];
	char id;
};
unsigned char Buffer[12][640];

extern struct projection *Proj;

int translate();
double get_dist();
double set_scale();
double rotate();
double strech();

main()
{
	int i,j,k;
	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;

	do
	{
		printf("\n   p -- set projection\n");
		printf("   l -- load constants\n");
		printf("   c -- check lat/lon to x/y and back again\n");
		printf("   g -- plot a graticule\n");
		printf("   P -- plot coastlines\n");
		printf("   d -- debug graticule\n");
		printf("   f -- fly\n");
		printf("   e -- expand last P\n");
		printf("   t -- translate xy to lat/lon\n");
		printf("\n   x -- exit program\n\n");

		ans=getch();

		if(ans=='p')
		{
			do
			{
				pnum=choose_projection(&Constant);
			}while(pnum<0);
/*			printf("\nProjection = %s\n\n",Proj[pnum].name);*/
			printf("Constant.proj = %c\n",Constant.proj);
		}
		if(ans=='t')
		{
			pnum=choose_projection(&Constant);
			if(pnum>=0)
			{
/*				printf("\nProjection = %s\n\n",Proj[pnum].name);*/
				printf("Constant.proj = %c\n",Constant.proj);
				printf("Give first guess at parameters.\n\n");
				get_constants(&Constant);
				translate(&Constant,Lat,Lon,X,Y,Tx,Ty,MAX_PTS);
			}
		}
		if(ans=='l')
		{
			get_constants(&Constant);
	}
		if(ans=='c')
		{
			printf("give lat and lon for check.\n");
			scanf("%lf%lf",&lat,&lon);
			ret=ll2xy(lat,lon,&x,&y,Constant);
			if(ret>=0)
			{
				printf("x = %lf  y = %lf\n",x,y);
				ret=xy2ll(x,y,&lat,&lon,Constant);
				if(ret>=0)
					printf("lat = %lf  lon = %lf\n",lat,lon);
				else
					printf("Error in function 'xy2ll()'.\n");
			}
			else
				printf("Error in function 'll2xy()'.\n");
		}
		if(ans=='g')
		{
			xc=0;yc=0;xs=640;ys=480;speed=8;
			printf("Give lat and lon of image center.\n");
			scanf("%lf%lf",&clat,&clon);
			limit_area(&xc,&yc,&xs,&ys,&speed,640,1,Buffer);
			horizon(xc,yc,xs,ys,clat,clon,Constant,D_BLUE,0);
			graticule(xc,yc,xs,ys,clat,clon,
				Constant,D_GREY,D_GREY,(float)10.0,0,(float)4.0,Constant);
			getch();
		}
		if(ans=='P')
		{
			xc=0;yc=0;xs=640;ys=480;speed=8;
			printf("Give lat and lon of image center.\n");
			scanf("%lf%lf",&clat,&clon);
			set_lut('c');
			paint_box(3,0,Buffer,0,0,640,480);
			limit_area(&xc,&yc,&xs,&ys,&speed,640,1,Buffer);
			horizon(xc,yc,xs,ys,clat,clon,Constant,D_BLUE,0);
			graticule(xc,yc,xs,ys,clat,clon,
				Constant,D_GREY,D_GREY,(float)10.0,0,(float)4.0,Constant);
			if(Constant.scale>=30.0)
				convert_bin(GREEN,"geo0",xc,yc,xs,ys,clat,clon,Constant);
			else
				convert_bin(GREEN,"geocoast",xc,yc,xs,ys,clat,clon,Constant);
			compress("project.seg",Buffer[0]);
			getch();
		}
		if(ans=='d')
		{
			xc=0;yc=0;xs=640;ys=480;speed=8;
			printf("Give lat and lon of image center.\n");
			scanf("%lf%lf",&clat,&clon);
			horizon(xc,yc,xs,ys,clat,clon,Constant,D_BLUE,1);
			ret=graticule(xc,yc,xs,ys,clat,clon,
				Constant,D_GREY,D_GREY,(float)10.0,1,(float)4.0,Constant);
			printf("graticule() = %d\n",ret);
			getch();
		}
		if(ans=='e')
		{
			set_lut('c');
			i=0;
			do
			{
				sprintf(filename,"proj%d.seg",i++);
				ret=expand(filename,Buffer[0]);
			}while(ret>=0);
			getch();
		}
		if(ans=='f')
		{
			i=0;
			fpin=fopen("project.ctl","rt");
			xc=yc=0;
			xs=640;ys=480;
			Constant.proj='P';
			while(fscanf(fpin,"%lf%lf%lf%lf%lf%lf%lf%lf",
				&Constant.scale,
				&Constant.H,
				&Constant.lat_1,
				&Constant.lon_0,
				&Constant.gamma,
				&Constant.omega,
				&clat,
				&clon)==8)
			{
/*				auto_set_constants(&Constant);*/
				paint_box(3,0,Buffer,0,0,640,480);
				horizon(xc,yc,xs,ys,clat,clon,Constant,D_BLUE,0);
				ret=graticule(xc,yc,xs,ys,clat,clon,
					Constant,D_GREY,D_GREY,(float)10.0,0,(float)4.0,Constant);
				if(ret<0)exit(0);
/*				ret=convert_bin(GREEN,"geocoast",xc,yc,xs,ys,clat,clon,Constant);
				if(ret<0)exit(0);
				sprintf(filename,"proj%d.seg",i++);
				compress(filename,Buffer[0]);*/
			}
			getch();
			fclose(fpin);
		}
	}while(ans!='x');
}

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

int video_on()

{
	int i,row;

	i=GetVideoBoardID();
	row=SetVideoMode(0x12);
	if(row==480)
	{
		ScreenXs=640;
		ScreenYs=480;
	}
	else if(row==400)
	{
		ScreenXs=640;
		ScreenYs=400;
	}
	else
	{
		row=SetVideoMode(0x10);
		if(row!=350)
		{
			SetVideoMode(0);
			printf("Could not boot color board.\n");
			exit(0);
		}
		ScreenXs=640;
		ScreenYs=350;
	}
	if(row<0)
	{
		SetVideoMode(0);
		printf("Could not boot board -- row = %d\n\n",row);
		exit(0);
	}
}

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

int video_off()

{
	SetVideoMode(0);
	return;
}





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

int translate(cc,lat,lon,x,y,tx,ty,max)

struct control *cc;
double *lat,*lon,*x,*y,*tx,*ty;
int max;

{
	int i,j,k;
	FILE *fpin,*fpout;
	char file[100],string[100];
	int num_pts=-1,val,ret;
	double dist,xb,yb,dnum,dscale=2.0,scale,d,d_up,d_down,error,dxcor=0.01;
	double theta=0.0,dtheta=.1,xcor=1.0;
	double llat,llon,lx,ly,mx,my,cost,sint;
	double xdist,ydist,dsc=R*PI/180.0;

	printf("proj = %c\n",cc->proj);
	printf("lon_0=%lf lon_1=%lf lon_2=%lf lon_p=%lf\n",
		cc->lon_0,cc->lon_1,cc->lon_2,cc->lon_p);
	printf("lat_0=%lf lat_1=%lf lat_2=%lf lat_p=%lf\n",
		cc->lat_0,cc->lat_1,cc->lat_2,cc->lat_p);
	printf("pixsize=%lf scale=%lf\n",cc->pixsize,cc->scale);
	printf("x0=%lf y0=%lf\n",cc->x0,cc->y0);
	printf("c=%lf C=%lf F=%lf H=%lf k=%lf\n",cc->c,cc->C,cc->F,cc->H,cc->k);
	printf("k0=%lf k1=%lf k2=%lf\n",cc->k0,cc->k1,cc->k2);
	printf("n=%lf P=%lf rho_0\n",cc->n,cc->P,cc->rho_0);
	printf("omega=%lf gamma=%lf horiz=%lf\n",cc->omega,cc->gamma,cc->horiz);
	getch();


	printf("Give input file name.  Type will be '.xyc'.\n");
	scanf("%s",file);
	strcpy(string,file);
	strcat(string,".xyc");
	fpin=fopen(string,"rt");
	if(!fpin)
	{
		printf("Could not open file '%s'.\n",string);
		return(-1);
	}
	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 %5.0lf %5.0lf\n",lat[i],lon[i],x[i],y[i]);
		i+=1;
	}
	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,6,xcor);
printf("Error after 'set_scale' = %lf\n",error);
	dscale=1.001;
	for(i=0;i<5;i++)
	{
		error=set_scale(num_pts,lat,lon,x,y,tx,ty,cc,&xb,&yb,dscale,theta,1,xcor);
		error=rotate(num_pts,lat,lon,x,y,tx,ty,cc,&xb,&yb,dtheta,&theta,1,xcor);
		error=strech(num_pts,lat,lon,x,y,tx,ty,cc,&xb,&yb,dxcor,theta,1,&xcor);
		dscale=(dscale-1.0)/10.0+1.0;
		dtheta/=10.0;
		dxcor/=10.0;
		printf("*");
	}
	printf("\n\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("  strech = %8.6lf\n",xcor);
	printf("rotation = %12.5lf deg.\n",theta/DEG_RAD);
	printf("tablet root mean error = %10.3lf mm\n",error/40.0);
	printf(" earth root mean error = %10.3lf km = %10.3lf mi\n",
		error*R/cc->scale,error*R/cc->scale*KM_MILE);
	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        ");
	printf("           errors (km)\n");
	printf("-- ------------------   -------------------  ");
	printf(" --------------------------\n");
	printf("      Lat      Lon         Lat       Lon     ");
	printf("    X        Y       Total\n");
	printf("-- -------- ---------   --------- ---------  ");
	printf(" -------  -------   -------\n");
	for(i=0;i<num_pts;i++)
	{
		lx=x[i];
		ly=y[i];
		lx=(lx+xb)*xcor;
		ly+=yb;
		mx=lx*cost-ly*sint;
		my=lx*sint+ly*cost;
		my=12200-my;
		llat=llon=0.0;
		ret=xy2ll(mx,my,&llat,&llon,*cc);
		if(ret>=0)
		{
			xdist=(llon-lon[i])*cos(lat[i])*dsc;
			ydist=(llat-lat[i])*dsc;
			dist=sqrt(xdist*xdist+ydist*ydist);
		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);
		}
	}

	strcpy(string,file);
	strcat(string,".raw");
	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;
			lx=(lx+xb)*xcor;
			ly+=yb;
			mx=lx*cost-ly*sint;
			my=lx*sint+ly*cost;
			my=12200-my;
			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);
	fclose(fpout);
}

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

int test_grat(num,lat,lon,x,y,cc)

int num;
double *lat,*lon,*x,*y;
struct control cc;

{
	int i,j,k;
	double scale=36.0;
	int ix,iy;

	for(i=0;i<num;i++)
	{
		ix=x[i]/scale;
		iy=y[i]/scale;
		cursor(3,ix,iy,255,7);
	}
}


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

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

int num_pts,times;
struct control *cc;
double *lat,*lon,*x,*y,*tx,*ty,*xbar,*ybar,dscale,*theta,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 strech(num_pts,lat,lon,x,y,tx,ty,cc,xbar,ybar,dscale,theta,
	times,xcor)

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

{
	int i,j,k;
	double dist,xb,yb,dnum,dxcor=dscale,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+dscale;
			d_up=get_dist(num_pts,lat,lon,x,y,tx,ty,cc,&xb,&yb,theta,*xcor);
			*xcor=xcor0+dscale;
			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(num_pts,lat,lon,x,y,tx,ty,cc,xbar,ybar,dscale,theta,
		times,xcor)

int num_pts,times;
struct control *cc;
double *lat,*lon,*x,*y,*tx,*ty,*xbar,*ybar,dscale,theta,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)/10.0+1.0;
	}
	*xbar=xb;
	*ybar=yb;
	return(d);
}

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

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

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

{
	int i,j,k;
	double xb=*xbar,yb=*ybar,dist,dnum=num_pts;
	double xx,yy;

	theta;

	xb=yb=0.0;
	for(i=0;i<num_pts;i++)
	{
		ll2xy(lat[i],lon[i],&xx,&yy,*cc);
		yy=12200-yy;
		tx[i]=xx*cos(theta)-yy*sin(theta);
		ty[i]=xx*sin(theta)+yy*cos(theta);
		xb+=(tx[i]-x[i])/dnum;
		yb+=(ty[i]-y[i])/dnum;
	}
	dist=0.0;
	for(i=0;i<num_pts;i++)
	{
		dist+=(tx[i]-(x[i]+xb)*xcor)*(tx[i]-(x[i]+xb)*xcor)/dnum+
				(ty[i]-(y[i]+yb))*(ty[i]-(y[i]+yb))/dnum;
	}
	*xbar=xb;
	*ybar=yb;
	return(sqrt(dist));
}