/* vmaputil.c

*/

/* ***** include ***** */

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

/* ***** define ***** */

#define NUMMAP 6
#define R      6367.40
#define PI		3.1415927
#define PI4		0.7853981
#define K0		0.9335351
#define DEG_RAD 0.01745329251				/* converts degrees to radians */
#define NUMCOAST 4
#define PAUSE 				printf("HIT ANY KEY TO CONTINUE\n");getch()

/* ***** variables ***** */

	/* ***** internal ***** */

	double MU_MinLon,MU_MaxLon;

	/* ***** external ***** */
	extern int Rxc,Ryc,Repeat;		/* multiple screens if Repeat>1 */

/* ***** arrays ***** */

	/* ***** internal ***** */

	char pstring[80];

	char *coastfile[10]={"coast.dat","namap.dat","samap.dat","eumap.dat",
		"aimap.dat","asmap.dat","afmap.dat","anmap.dat",
		"country.dat","eus.dat"};

	/* ***** external ***** */

	extern int Font[128][25];

/* ***** stuctures ***** */

	/* ***** external ***** */

	extern struct map_screen
	{
		int xc,yc,xs,ys;
	}map_pos;

	/* ***** definitions ***** */

	/* ***** declarations ***** */

	struct control
	{
		char type;				/* s=polar stereo  p=plate carree  m=Mercator */
		int  tnum;				/* projection number */
		char *name;				/* projection name */
		double clat,clon; 	/* image center */
		double pixsize;		/* reference pixel size at some point */
		double scale,scale2;	/* scale factor (2 is for plate carree) */
		double lon0; 			/* prime longitude (0 is down or up for pstereo */
		double x0,y0;			/* x,y of upper left corner */
		double F,n,rho0;		/* Lambert conformal constants */
		double C,lat1,lat2;  /* Albers const. and stand. parallels */
	}cc;

	struct maps
	{
		char name[40];			/* projection name */
		char id;					/* projection identifier */
		char pixsize[40];		/* location of pixel size */
	}mm[NUMMAP]=
		{
			{"polar stereographic",'s',"at pole"},
			{"modified plate carree",'p',"at image center"},
			{"Mercator",'m',"at equator"},
			{"Lambert Conformal",'l',"at image center"},
			{"Albers Equal-Area",'a',"at image center"},
			{"sinusoidal",'S',"at image center"}
		};

/* ***** functions ***** */

	/* ***** external ***** */

	void plot_font_ovl();
	void cplotln();
	double tan(),exp(),pow();

	/* ***** internal ***** */

	void 	auto_set_cc( char , double * );
	void 	cplotln(int,int,int,int,int,int);
	int 	d_off_image( double , double );
	int 	d_on_image( double , double );
	int 	i_off_image( int , int );
	int 	i_on_image( int , int );
	int 	image_max_min( double * , double * , double * , double * );
	void	 lat_lon_of_box( void );
	void	 ll_xy( double * , double * , double , double );
	void	 plot_coast( int , int , int , int );
	void	 plot_graticule( int * , int , int , double , int);
	void	 plot_graticule_frac( int * , int , int , double );
	void	 print_cc( FILE * , int );
	void	 set_cc( void );
	void	 xy_ll( double , double , double * , double * );



/********************************************************************
**
**	sets all the constants needed to map data
**
**	You can call it instead of set_cc()
**
**		    type -- defines projection
**		indat[0] -- pixsize -- projection of interest
**		indat[1] -- clat 	  -- central latitude
**		indat[2] -- clon	  -- central longitude
**		indat[3] -- lon0    -- vertical graticule 'a','S' & 's' only
**		indat[4] -- lat1    -- 1st standard parallel 'a' only
**		indat[5] -- lat2    -- 2nd standard parallel 'a' only
**
********************************************************************* */

void auto_set_cc(type,indat)

char type;
double *indat;

{
	int i;
	int valid=0;
	double x,y,tval,sign=1.0;
	double mlat1,mlat2,mlon;
	double xc=map_pos.xc,yc=map_pos.yc;
	double xm=map_pos.xc+map_pos.xs/2,ym=map_pos.yc+map_pos.ys/2;

	do
	{
		for(i=0;i<NUMMAP;i++)
			if(type==mm[i].id)
			{
				valid=1;
				cc.type=type;
				cc.tnum=i;
				cc.name=mm[i].name;
			}
	}while(valid==0&&type!='x');
	if(type=='x')return;

	cc.pixsize=indat[0];
	cc.clat=indat[1];
	cc.clon=indat[2];

	if(cc.type=='S')
	{
		cc.lon0=indat[3];
		MU_MaxLon=cc.clon+180;
		MU_MinLon=cc.clon-180;
		cc.scale=R*PI/180.0/cc.pixsize;	/* pixels/deg of lat */
		cc.x0=0.0;
		cc.y0=0.0;
		ll_xy(&x,&y,cc.clat,cc.clon);
		cc.x0=x-xm;
		cc.y0=y-xm;
		xy_ll((double)256.5,(double)256.5,&mlat1,&mlon);
		xy_ll((double)256.5,(double)255.5,&mlat2,&mlon);
	}
	if(cc.type=='a')
	{
		cc.lat1=indat[4];
		cc.lat2=indat[5];
		cc.lon0=indat[3];
		MU_MaxLon=cc.clon+180.0;
		MU_MinLon=cc.clon-180.0;
		cc.scale=R/cc.pixsize;
		if(cc.clat<0.0)sign=-1.0;
		cc.lat1=cc.lat1*sign;
		cc.lat2=cc.lat2*sign;
		cc.n=(sin(cc.lat1*DEG_RAD)+sin(cc.lat2*DEG_RAD))/2.0;
		cc.C=cos(cc.lat1*DEG_RAD)*cos(cc.lat1*DEG_RAD)+
				2.0*cc.n*sin(cc.lat1*DEG_RAD);
		cc.rho0=R*sqrt(cc.C-2.0*cc.n*sin(cc.clat*DEG_RAD))/cc.n;
		cc.x0=xc;
		cc.y0=yc;
		ll_xy(&x,&y,cc.clat,cc.clon);
		cc.x0=x-xm;
		cc.y0=y-ym;
	}
	if(cc.type=='s')	/* polar stereographic */
	{
		MU_MaxLon=180.0;
		MU_MinLon=-180.0;
		cc.scale=2.0*R*K0/cc.pixsize;
		if(cc.clat>0.0)
			cc.lat1=90.0;
		else
			cc.lat1=-90.0;
		cc.lon0=indat[3];
		cc.x0=0;
		cc.y0=0;
		ll_xy(&x,&y,cc.clat,cc.clon);
		cc.x0=x-xm;
		cc.y0=y-ym;
		xy_ll((double)256.5,(double)256.5,&mlat1,&mlon);
		xy_ll((double)256.5,(double)255.5,&mlat2,&mlon);
	}
	if(cc.type=='p')	/* plate carree */
	{
		MU_MaxLon=cc.clon+180.0;
		MU_MinLon=cc.clon-180.0;
		cc.scale=R*PI/180.0/cc.pixsize;	/* pixels/deg of lat */
		cc.scale2=cc.scale*cos(cc.clat*DEG_RAD); /* pixels/deg of lon */
		cc.lon0=cc.clon;
		cc.lat1=cc.clat;
		cc.x0=-xm;
		cc.y0=-ym;
	}
	if(cc.type=='m')	/* Mercator */
	{
		MU_MaxLon=cc.clon+180.0;
		MU_MinLon=cc.clon-180.0;
		cc.scale=R*PI/180.0/cc.pixsize;	/* pixels/deg of lat */
		cc.lon0=0;
		cc.lat1=0;
		cc.x0=0;
		cc.y0=0;
		ll_xy(&x,&y,cc.clat,cc.clon);
		cc.x0=x-xm;
		cc.y0=y-ym;
	}	
	if(cc.type=='l')
	{
		MU_MaxLon=cc.clon+180.0;
		MU_MinLon=cc.clon-180.0;
		if(cc.clat<0.0)sign=-1.0;
		cc.scale=R/cc.pixsize;
		cc.lon0=cc.clon;
		cc.lat1=cc.clat*sign;
		cc.n=sin(cc.lat1*DEG_RAD);
		tval=pow(tan(PI4+cc.lat1/2.0*DEG_RAD),cc.n);
		cc.F=cos(cc.lat1*DEG_RAD)*tval/cc.n;
		cc.rho0=R*cc.F/tval;
		cc.x0=0;
		cc.y0=0;
		ll_xy(&x,&y,cc.clat,cc.clon);
		cc.x0=x-xm;
		cc.y0=y-ym;
	}
}



/*********************************************************************
**
**	checks if double x,y pixel is on image:
**
**	returns:  'false' -- 0 -- off image
**				 'true'  -- 1 -- on image
**
********************************************************************** */

int d_on_image(x,y)

double x,y;

{
	double x1=map_pos.xc,y1=map_pos.yc;
	double x2=map_pos.xc+map_pos.xs-1,y2=map_pos.yc+map_pos.ys-1;

	if(x<x1||x>x2||y<y1||y>y2)
		return(0);
	else
		return(1);
}

/*********************************************************************
**
**	checks if double x,y pixel is off image:
**
**	returns:  'false' -- 0 -- on image
**				 'true'  -- 1 -- off image
**
********************************************************************** */

int d_off_image(x,y)

double x,y;

{
	double x1=map_pos.xc,y1=map_pos.yc;
	double x2=map_pos.xc+map_pos.xs-1,y2=map_pos.yc+map_pos.ys-1;

	if(x<x1||x>x2||y<y1||y>y2)
		return(1);
	else
		return(0);
}


/*********************************************************************
**
**	checks if int x,y pixel is on image:
**
**	returns:  'false' -- 0 -- off image
**				 'true'  -- 1 -- on image
**
********************************************************************** */

int i_on_image(x,y)

int x,y;

{
	int x1=map_pos.xc,y1=map_pos.yc;
	int x2=map_pos.xc+map_pos.xs-1,y2=map_pos.yc+map_pos.ys-1;

	if(x<x1||x>x2||y<y1||y>y2)
		return(0);
	else
		return(1);
}

/*********************************************************************
**
**	checks if int x,y pixel is off image:
**
**	returns:  'false' -- 0 -- on image
**				 'true'  -- 1 -- off image
**
********************************************************************** */

int i_off_image(x,y)

int x,y;

{
	int x1=map_pos.xc,y1=map_pos.yc;
	int x2=map_pos.xc+map_pos.xs-1,y2=map_pos.yc+map_pos.ys-1;

	if(x<x1||x>x2||y<y1||y>y2)
		return(1);
	else
		return(0);
}


/********************************************************************
**
**	plots lat/lon on overlay and lables them
**
********************************************************************* */

void plot_graticule_frac(font,lbrite,nbrite,del)

int *font,lbrite,nbrite;
double del;

{
	int l,ll;
	char string[80];
	double latmax,latmin,lonmax,lonmin;
	double x,y,lat,lon,plon;
	double x1=map_pos.xc+15.0,y1=map_pos.yc+10.0;
	double x2=map_pos.xc+map_pos.xs-16.0,y2=map_pos.yc+map_pos.ys-11.0;
	int ix1,iy1,ix2,iy2,l1,l2,ld,v;
	int pole;
	double dll;
	char *format;

	format="%0.1lf";
	if(del<0.3)
		format="%0.2lf";
	if(del<0.03)
		format="%0.3lf";
	if(del<0.003)
		format="%0.4lf";
	if(del<0.0003)
		format="%0.5lf";

	pole=image_max_min(&latmax,&latmin,&lonmax,&lonmin);

	ix1=lonmin/del;
	lon=ix1;
	lon*=del;
	ll=(latmax+latmin)/2/del;
	dll=ll;
	dll*=del;
	dll+=del/2;

	if(cc.type=='S')
	{
		do		/* plot longitude */
		{
			lat=latmax;
			do
			{
				ll_xy(&x,&y,lat,lon);
				ix1=x;iy1=y;
				ll_xy(&x,&y,lat-del/4,lon);
				ix2=x;iy2=y;
				cplotln(3,ix1,iy1,ix2,iy2,nbrite);
				lat-=del/4;
			}while(lat>=latmin+del/4);

			ll_xy(&x,&y,dll,lon);
			if(x>x1&&x<x2&&y>y1&&y<y2)
			{
				plon=lon;
				if(plon>180)plon-=360.0;
				if(plon<-180)plon+=360.0;
				sprintf(string,format,plon);
				ix1=x-strlen(string)*5/2;
				iy1=y;
				if(nbrite>0)
					plot_font_ovl(nbrite,ix1,iy1,string,3,font);
			}
			lon+=del;
		}while(lon<=lonmax);	

	}
	else
		do		/* plot longitude */
		{
			ll_xy(&x,&y,latmax,lon);
			ix1=x;iy1=y;
			ll_xy(&x,&y,latmin,lon);
			ix2=x;iy2=y;
			cplotln(3,ix1,iy1,ix2,iy2,nbrite);
			ll_xy(&x,&y,dll,lon);
			if(x>x1&&x<x2&&y>y1&&y<y2)
			{
				plon=lon;
				if(plon>180.0)plon-=360.0;
				if(plon<-180.0)plon+=360.0;
				sprintf(string,format,plon);
				ix1=x-strlen(string)*5/2;
				iy1=y;
				if(nbrite>0)
					plot_font_ovl(nbrite,ix1,iy1,string,3,font);
			}
			lon+=del;
		}while(lon<=lonmax);	

	ll=(lonmax+lonmin)/2/del;
	dll=ll;
	dll*=del;
	dll+=del/2;
	if(cc.type=='m'||cc.type=='p')	/* straight lines of latitude */
	{
		ix1=latmin/del;
		lat=ix1;
		lat*=del;
		do
		{
			ll_xy(&x,&y,lat,lonmax);
			iy1=y;
			cplotln(3,map_pos.xc,iy1,map_pos.xc+map_pos.xs-1,iy1,nbrite);
			ll_xy(&x,&y,lat,dll);
			if(x>x1&&x<x2&&y>y1&&y<y2)
			{
				sprintf(string,format,lat);
				ix1=x;
				iy1=y-5;
				if(nbrite>0)
					plot_font_ovl(nbrite,ix1,iy1,string,3,font);
			}
			lat+=del;
		}while(lat<=latmax);	
	}
	else										/* curved lines of latitude */
	{
		if(pole==1)
		{
			l1=MU_MinLon;
			l2=MU_MaxLon;
			l2-=1;
		}
		else
		{
			ld=1.0/del;
			l1=lonmin/del-1;
			l2=lonmax/del;
			l2+=1;
		}
		ix1=latmin/del-1;
		lat=ix1;
		lat*=del;
		do
		{
			for(l=l1;l<l2;l++)
			{
				lon=l;
				lon*=del;
				ll_xy(&x,&y,lat,lon);
				ix1=x;iy1=y;
				lon+=del;
				ll_xy(&x,&y,lat,lon);
				ix2=x;iy2=y;
				cplotln(3,ix1,iy1,ix2,iy2,nbrite);
			}
			ll_xy(&x,&y,lat,dll);
			if(x>x1&&x<x2&&y>y1&&y<y2)
			{
				sprintf(string,format,lat);
				ix1=x-strlen(string)*5/2;
				iy1=y;
				if(nbrite>0)
					plot_font_ovl(nbrite,ix1,iy1,string,3,font);
			}
			lat+=del;
		}while(lat<=latmax);	
	}
	plot_graticule(font,lbrite,0,(double)1.0,0);
}

/********************************************************************
**
**	plots lat/lon on overlay and lables them
**
********************************************************************* */

void plot_graticule(font,lbrite,nbrite,del,iop)

int *font,lbrite,nbrite,iop;
double del;

{
	int l,ll;
	char string[80];
	double latmax,latmin,lonmax,lonmin;
	double x,y,lat,lon,plon;
	double x1=map_pos.xc+15.0,y1=map_pos.yc+10.0;
	double x2=map_pos.xc+map_pos.xs-16.0,y2=map_pos.yc+map_pos.ys-11.0;
	int ix1,iy1,ix2,iy2,l1,l2,ld,v;
	int pole;
	FILE *fp;

	if(iop==1)
		fp=fopen("graticul.out","wt");

	if(del<1.0)
	{
		plot_graticule_frac(font,lbrite,nbrite,del);
		return;
	}

	pole=image_max_min(&latmax,&latmin,&lonmax,&lonmin);

	ix1=lonmin/del;
	lon=ix1;
	lon*=del;
	ll=(latmax+latmin)/2/del;
	ll*=del;
	ll+=del/2;

	if(cc.type=='S')
	{
		do		/* plot longitude */
		{
			lat=latmax;
			do
			{
				ll_xy(&x,&y,lat,lon);
				ix1=x;iy1=y;
				ll_xy(&x,&y,lat-del/4,lon);
				ix2=x;iy2=y;
				cplotln(3,ix1,iy1,ix2,iy2,lbrite);
				lat-=del/4;
			}while(lat>=latmin+del/4);

			ll_xy(&x,&y,(double)ll,lon);
			if(x>x1&&x<x2&&y>y1&&y<y2)
			{
				plon=lon;
				if(plon>180)plon-=360.0;
				if(plon<-180)plon+=360.0;
				sprintf(string,"%4.0lf",plon);
				ix1=x-strlen(string)*5/2;
				iy1=y;
				if(nbrite>0)
					plot_font_ovl(nbrite,ix1,iy1,string,3,font);
			}
			lon+=del;
		}while(lon<=lonmax);	

	}
	else
		do		/* plot longitude */
		{
			ll_xy(&x,&y,latmax,lon);
			ix1=x;iy1=y;
			ll_xy(&x,&y,latmin,lon);
			ix2=x;iy2=y;
			cplotln(3,ix1,iy1,ix2,iy2,lbrite);
			ll_xy(&x,&y,(double)ll,lon);
			if(x>x1&&x<x2&&y>y1&&y<y2)
			{
				plon=lon;
				if(plon>180.0)plon-=360.0;
				if(plon<-180.0)plon+=360.0;
				sprintf(string,"%4.0lf",plon);
				ix1=x-strlen(string)*5/2;
				iy1=y;
				if(nbrite>0)
					plot_font_ovl(nbrite,ix1,iy1,string,3,font);
				if(iop==1)
				{
					fprintf(fp,"%s  %d %d\n",string,ix1,iy1);
					fprintf(fp,"0\n%4d %4d  %4d %4d//\n",ix1,iy1,ix2,iy2);
				}
			}
			lon+=del;
		}while(lon<=lonmax);	

	ll=(lonmax+lonmin)/2/del;
	ll*=del;
	ll+=del/2;
	if(cc.type=='m'||cc.type=='p')	/* straight lines of latitude */
	{
		ix1=latmin/del;
		lat=ix1;
		lat*=del;
		do
		{
			ll_xy(&x,&y,lat,lonmax);
			iy1=y;
			cplotln(3,map_pos.xc,iy1,map_pos.xc+map_pos.xs-1,iy1,lbrite);
			ll_xy(&x,&y,lat,(double)ll);
			if(x>x1&&x<x2&&y>y1&&y<y2)
			{
				sprintf(string,"%4.0lf",lat);
				ix1=x;
				iy1=y-5;
				if(nbrite>0)
					plot_font_ovl(nbrite,ix1,iy1,string,3,font);
				if(iop==1)
				{
					fprintf(fp,"%s  %d %d\n",string,ix1,iy1);
					fprintf(fp,"1\n%4d %4d  %4d %4d//\n",
							map_pos.xc,iy1,map_pos.xc+map_pos.xs-1,iy1);
				}
			}
			lat+=del;
		}while(lat<=latmax);	
	}
	else										/* curved lines of latitude */
	{
		if(pole==1)
		{
			l1=MU_MinLon;
			l2=MU_MaxLon;
			l2-=1;
		}
		else
		{
			ld=del;
			l1=lonmin/del-1;
			l1*=ld;
			l2=lonmax/del;
			l2+=1;
			l2*=ld;
		}
		ix1=latmin/del;
		lat=ix1;
		lat*=del;
		do
		{
			for(l=l1;l<l2;l++)
			{
				lon=l;
				ll_xy(&x,&y,lat,lon);
				ix1=x;iy1=y;
				lon=l+1;
				ll_xy(&x,&y,lat,lon);
				ix2=x;iy2=y;
				cplotln(3,ix1,iy1,ix2,iy2,lbrite);
			}
			ll_xy(&x,&y,lat,(double)ll);
			if(x>x1&&x<x2&&y>y1&&y<y2)
			{
				sprintf(string,"%4.0lf",lat);
				ix1=x-strlen(string)*5/2;
				iy1=y;
				if(nbrite>0)
					plot_font_ovl(nbrite,ix1,iy1,string,3,font);
			}
			lat+=del;
		}while(lat<=latmax);	
	}
}


/*********************************************************************
**
**	finds max/min for plotting graticule in function above
**
********************************************************************** */

int image_max_min(latx,latn,lonx,lonn)

double *latx,*latn,*lonx,*lonn;	/* pointers to max/min lat/lon */

{
	double latmax,latmin,lonmax,lonmin;
	double x,y,lat,lon;
	double x1=map_pos.xc,y1=map_pos.yc;
	double x2=map_pos.xc+map_pos.xs-1,y2=map_pos.yc+map_pos.ys-1;
	double xm=map_pos.xc+map_pos.xs/2.0,ym=map_pos.yc+map_pos.ys/2.0;
	int pole=0;

	x=x1;y=y1;
	xy_ll(x,y,&lat,&lon);
	latmax=latmin=lat;
	lonmax=lonmin=lon;
	x=xm;y=y1;
	xy_ll(x,y,&lat,&lon);
	if(lat>latmax)latmax=lat;if(lat<latmin)latmin=lat;
	if(lon>lonmax)lonmax=lon;if(lon<lonmin)lonmin=lon;
	x=x2;y=y1;
	xy_ll(x,y,&lat,&lon);
	if(lat>latmax)latmax=lat;if(lat<latmin)latmin=lat;
	if(lon>lonmax)lonmax=lon;if(lon<lonmin)lonmin=lon;
	x=x1;y=ym;
	xy_ll(x,y,&lat,&lon);
	if(lat>latmax)latmax=lat;if(lat<latmin)latmin=lat;
	if(lon>lonmax)lonmax=lon;if(lon<lonmin)lonmin=lon;
	x=x2;y=ym;
	xy_ll(x,y,&lat,&lon);
	if(lat>latmax)latmax=lat;if(lat<latmin)latmin=lat;
	if(lon>lonmax)lonmax=lon;if(lon<lonmin)lonmin=lon;
	x=x1;y=y2;
	xy_ll(x,y,&lat,&lon);
	if(lat>latmax)latmax=lat;if(lat<latmin)latmin=lat;
	if(lon>lonmax)lonmax=lon;if(lon<lonmin)lonmin=lon;
	x=xm;y=y2;
	xy_ll(x,y,&lat,&lon);
	if(lat>latmax)latmax=lat;if(lat<latmin)latmin=lat;
	if(lon>lonmax)lonmax=lon;if(lon<lonmin)lonmin=lon;
	x=x2;y=y2;
	xy_ll(x,y,&lat,&lon);
	if(lat>latmax)latmax=lat;if(lat<latmin)latmin=lat;
	if(lon>lonmax)lonmax=lon;if(lon<lonmin)lonmin=lon;

	if(latmax>90.0)latmax=90.0;
	if(latmin<-90.0)latmin=-90.0;

	if(cc.type=='l'||cc.type=='s'||cc.type=='a'||cc.type=='S')/*pole in image*/
	{
		if(cc.clat>0.0)	/* N hemisphere */
		{
			lat=90.0;lon=0;
			ll_xy(&x,&y,lat,lon);
			if(d_on_image(x,y))
			{
				pole=1;
				latmax=90.0;
				lonmin=MU_MinLon;
				lonmax=MU_MaxLon;
			}
		}
		else				/* S hemisphere */
		{
			lat=-90.0;lon=0;
			ll_xy(&x,&y,lat,lon);
			if(d_on_image(x,y))
			{
				pole=1;
				latmin=-90.0;
				lonmin=MU_MinLon;
				lonmax=MU_MaxLon;
			}
		}
	}
	if(cc.type=='S')
	{
		lonmin=MU_MinLon;
		lonmax=MU_MaxLon;
	}

	*latx=latmax;
	*latn=latmin;
	*lonx=lonmax;
	*lonn=lonmin;

	return(pole);
}



/********************************************************************
**
**	finds lat/lon of image corners, middle, and mid-sides
**
********************************************************************* */

void lat_lon_of_box()

{
	double lat,lon,x,y;
	double x1=map_pos.xc,y1=map_pos.yc;
	double x2=map_pos.xc+map_pos.xs-1,y2=map_pos.yc+map_pos.ys-1;
	double xm=map_pos.xc+map_pos.xs/2.0,ym=map_pos.yc+map_pos.ys/2.0;

	if(cc.type=='\0')
		set_cc();

	x=x1;y=y1;
	xy_ll(x,y,&lat,&lon);
	printf("\n\n%6.2lf,%7.2lf",lat,lon);
	x=xm;y=y1;
	xy_ll(x,y,&lat,&lon);
	printf("     %6.2lf,%7.2lf",lat,lon);
	x=x2;y=y1;
	xy_ll(x,y,&lat,&lon);
	printf("     %6.2lf,%7.2lf\n",lat,lon);
	x=x1;y=ym;
	xy_ll(x,y,&lat,&lon);
	printf("%6.2lf,%7.2lf",lat,lon);
	x=xm;y=ym;
	xy_ll(x,y,&lat,&lon);
	printf("     %6.2lf,%7.2lf",lat,lon);
	x=x2;y=ym;
	xy_ll(x,y,&lat,&lon);
	printf("     %6.2lf,%7.2lf\n",lat,lon);
	x=x1;y=y2;
	xy_ll(x,y,&lat,&lon);
	printf("%6.2lf,%7.2lf",lat,lon);
	x=xm;y=y2;
	xy_ll(x,y,&lat,&lon);
	printf("     %6.2lf,%7.2lf",lat,lon);
	x=x2;y=y2;
	xy_ll(x,y,&lat,&lon);
	printf("     %6.2lf,%7.2lf\n\n",lat,lon);
	PAUSE;
}

/********************************************************************
**
**	plots coast lines: 
**
**  0 = coastlines   1 = country boundaries
**
**		iop = 1 puts points on map into submap.dat
**
********************************************************************* */

void plot_coast(filenum,brite,iop,maxdist)

int filenum,brite,iop,maxdist;

{
	int i,j,k;
	FILE *fp,*fpout;
	double lat,lon,x,y;
	long ix,iy,ixo=0,iyo=0;
	int ll[2],val;

	maxdist*=maxdist;

	fp=fopen(coastfile[filenum],"rb");
	if(iop==1)fpout=fopen("submap.dat","wb");
	while(fread((char *)ll,sizeof(int),2,fp)>1)
	{
		lat=ll[0];
		lon=ll[1];
		lat/=100.0;
		lon/=100.0;
		if(lon>180.0)lon-=360.0;
		ll_xy(&x,&y,lat,lon);
		if(d_on_image(x,y))
		{
			ix=x;iy=y;
			if((ix-ixo)*(ix-ixo)<maxdist&&(iy-iyo)*(iy-iyo)<maxdist)
				cplotln(3,(int)ixo,(int)iyo,(int)ix,(int)iy,brite);
			ixo=ix;iyo=iy;
			if(iop==1)
				fwrite((char *)ll,sizeof(int),2,fpout);
		}
	}
	fclose(fp);
	if(iop==1)fclose(fpout);	
}



/********************************************************************
**
**	converts lat/lon to image x,y 
**
**	calls set_cc if needed
**
********************************************************************* */

void ll_xy(x,y,lat,lon)

double *x,*y,lat,lon;

{
	double rho,c,theta,tval,sign=1.0,scale=Repeat;

	if(cc.type=='\0')
		set_cc();

	if(lon>MU_MaxLon)lon-=360.0;
	if(lon<MU_MinLon)lon+=360.0;	

	if(cc.type=='s')
	{
		if(cc.lat1==90.0)
		{
			*x=cc.scale*tan(PI4-lat*DEG_RAD/2.0)*sin((lon-cc.lon0)*DEG_RAD);
			*y=cc.scale*tan(PI4-lat*DEG_RAD/2.0)*cos((lon-cc.lon0)*DEG_RAD);
		}
		else
		{
			*x=cc.scale*tan(PI4+lat*DEG_RAD/2.0)*sin((lon-cc.lon0)*DEG_RAD);
			*y=-cc.scale*tan(PI4+lat*DEG_RAD/2.0)*cos((lon-cc.lon0)*DEG_RAD);
		}

		*x-=cc.x0;
		*y-=cc.y0;
	}
	else if(cc.type=='p')
	{
		*x=(lon-cc.lon0)*cc.scale2-cc.x0;
		*y=(cc.lat1-lat)*cc.scale-cc.y0;
	}
	else if(cc.type=='m')
	{
		*x=(lon-cc.lon0)*cc.scale-cc.x0;
		*y=-cc.scale*log(tan(PI4+lat*DEG_RAD/2.0))/DEG_RAD-cc.y0;
	}
	else if(cc.type=='l')
	{
		if(cc.clat<0.0)sign=-1.0;
		lat*=sign;
		tval=pow(tan(PI4+lat/2.0*DEG_RAD),cc.n);
		rho=cc.scale*cc.F/tval;
		theta=cc.n*(lon-cc.lon0)*DEG_RAD;
		*x=rho*sin(theta)-cc.x0;
		*y=(cc.rho0-rho*cos(theta))*(-1.0)*sign-cc.y0;
	}
	else if(cc.type=='a')
	{
		if(cc.clat<0.0)sign=-1.0;
		lat*=sign;
		rho=cc.scale*sqrt(cc.C-2.0*cc.n*sin(lat*DEG_RAD))/cc.n;
		theta=cc.n*(lon-cc.lon0)*DEG_RAD;
		*x=rho*sin(theta)-cc.x0;
		*y=(cc.rho0-rho*cos(theta))*(-1.0)*sign-cc.y0;
	}
	else if(cc.type=='S')
	{
		*x=cc.scale*(lon-cc.clon)*cos(lat*DEG_RAD)-cc.x0;
		*y=(cc.clat-lat)*cc.scale-cc.y0;
	}
	if(Repeat>1)
	{
		*x=*x*scale-(double)Rxc;
		*y=*y*scale-(double)Ryc;
	}
}




/********************************************************************
**
**	converts image x,y to lat/lon
**
**	calls set_cc if needed
**
********************************************************************* */

void xy_ll(x,y,lat,lon)

double x,y,*lat,*lon;

{
	double rho,c,theta,sign=1.0,xx,yy,scale=Repeat;

	if(Repeat>1)
	{
		x=(x+(double)Rxc)/scale;
		y=(y+(double)Ryc)/scale;
	}
	if(cc.type=='\0')
		set_cc();

	if(cc.type=='s')
	{
		if(cc.clat<0.0)sign=-1.0;
		x+=cc.x0;
		y+=cc.y0;
		rho=sqrt(x*x+y*y)/cc.scale;
		c=2.0*atan(rho);
		*lat=asin(cos(c))/DEG_RAD*sign;
		if(cc.clat>0.0)
			*lon=cc.lon0+atan2(x,-y)/DEG_RAD;
		else
			*lon=cc.lon0+atan2(x,-y)/DEG_RAD;
	}
	else if(cc.type=='p')
	{
		*lon=(x+cc.x0)/cc.scale2+cc.lon0;
		*lat=-(y+cc.y0)/cc.scale+cc.lat1;
	}
	else if(cc.type=='m')
	{
		*lon=(x+cc.x0)/cc.scale+cc.lon0;
		*lat=90.0-2.0*atan(exp((y+cc.y0)/cc.scale*DEG_RAD))/DEG_RAD;
	}
	else if(cc.type=='l')
	{
		if(cc.clat<0.0)sign=-1.0;
		xx=x+cc.x0;
		yy=-(y+cc.y0);
		yy*=sign;
		rho=sqrt(xx*xx+(cc.rho0-yy)*(cc.rho0-yy));
		theta=atan2(xx,cc.rho0-yy);
		*lat=2.0*atan(pow(cc.scale*cc.F/rho,1.0/cc.n))-PI4*2.0;
		*lat/=DEG_RAD*sign;
		*lon=theta/cc.n/DEG_RAD+cc.lon0;
	}
	else if(cc.type=='a')
	{
		if(cc.clat<0.0)sign=-1.0;
		xx=x+cc.x0;
		yy=-(y+cc.y0);
		yy*=sign;
		rho=sqrt(xx*xx+(cc.rho0-yy)*(cc.rho0-yy));
		theta=atan2(xx,cc.rho0-yy);
		*lat=asin((cc.C-(rho*cc.n/cc.scale)*(rho*cc.n/cc.scale))/(2.0*cc.n));
		*lat/=DEG_RAD*sign;
		*lon=theta/cc.n/DEG_RAD+cc.lon0;
	}
	else if(cc.type=='S')
	{
		*lat=-(y+cc.y0)/cc.scale+cc.clat;
		*lon=cc.lon0+(x+cc.x0)/(cc.scale*cos(*lat*DEG_RAD));
	}
	if(*lon>MU_MaxLon)*lon-=360.0;
	if(*lon<MU_MinLon)*lon+=360.0;	
}




/********************************************************************
**
**	prints struct cc.xxx	on screen and into a file if option = 1
**
**	if option != 1 or *fpout = null it only prints on screen
**
********************************************************************* */

void print_cc(fpout,iop)

FILE *fpout;
int iop;

{
	if(cc.type=='\0')
		set_cc();

	printf("\n\nMap data -- struct control cc\n\n");
	printf("   cc.type = %c\n",cc.type);
	printf("   cc.tnum = %d\n",cc.tnum);
	printf("   cc.name = %s\n",cc.name);
	printf("   cc.clat = %lf    cc.clon = %lf\n",cc.clat,cc.clon);
	printf("   cc.pixsize = %lf\n",cc.pixsize);
	printf("   cc.scale = %lf   cc.scale2 = %lf\n",cc.scale,cc.scale2);
	printf("   cc.lon0 = %lf    cc.lat1 = %lf\n",cc.lon0,cc.lat1);
	printf("   cc.x0   = %lf    cc.y0   = %lf\n",cc.x0,cc.y0);
	if(cc.type=='l')
		printf("   cc.n = %lf   cc.F = %lf   cc.rho0 = %lf\n",
			cc.n,cc.F,cc.rho0);
	if(cc.type=='a')
		printf("   cc.n = %lf   cc.C = %lf   cc.rho0 = %lf   lat2 = %lf\n",
			cc.n,cc.C,cc.rho0,cc.lat2);

	if(iop==1&&fpout)
	{
		printf("\n");
		fprintf(fpout,"\n");
		fprintf(fpout,"\n\nMap data -- struct control cc\n\n");
		fprintf(fpout,"   cc.type = %c\n",cc.type);
		fprintf(fpout,"   cc.tnum = %d\n",cc.tnum);
		fprintf(fpout,"   cc.name = %s\n",cc.name);
		fprintf(fpout,"   cc.clat = %lf    cc.clon = %lf\n",cc.clat,cc.clon);
		fprintf(fpout,"   cc.pixsize = %lf\n",cc.pixsize);
		fprintf(fpout,"   cc.scale = %lf   cc.scale2 = %lf\n",cc.scale,cc.scale2);
		fprintf(fpout,"   cc.lon0 = %lf    cc.lat1 = %lf\n",cc.lon0,cc.lat1);
		fprintf(fpout,"   cc.x0   = %lf    cc.y0   = %lf\n",cc.x0,cc.y0);
		if(cc.type=='l')
			fprintf(fpout,"   cc.n = %lf   cc.F = %lf   cc.rho0 = %lf\n",
				cc.n,cc.F,cc.rho0);
	}

	PAUSE;
}




/********************************************************************
**
**	sets all the constants needed to map data
**
**	You can call it or it will be called by ll_xy or xy_ll if needed
**
********************************************************************* */

void set_cc()

{
	int i;
	char type;
	int valid=0;
	double x,y,tval,sign=1.0;
	double mlat1,mlat2,mlon;
	double xc=map_pos.xc+15.0,yc=map_pos.yc+10.0;
	double xm=map_pos.xc+map_pos.xs/2.0,ym=map_pos.yc+map_pos.ys/2.0;

	printf("\n\nEnter projection:\n\n");
	for(i=0;i<NUMMAP;i++)
		printf("   %c -- %s\n",mm[i].id,mm[i].name);
	printf("\n   x -- exit without setting map control\n");

	do
	{
		type=getch();
		for(i=0;i<NUMMAP;i++)
			if(type==mm[i].id)
			{
				valid=1;
				cc.type=type;
				cc.tnum=i;
				cc.name=mm[i].name;
			}
	}while(valid==0&&type!='x');
	if(type=='x')return;

	printf("\nMap type = %s\n\n",cc.name);

	do
	{
		printf("\nGive pixel size %s in km.\n",mm[cc.tnum].pixsize);
		scanf("%lf",&cc.pixsize);
		if(cc.pixsize<=0.0)
			printf("\n\nPIXEL SIZE INVALID !  TRY AGAIN !\n\n");
	}while(cc.pixsize<=0.0);

	printf("\nPixel size %s = %lf km.\n",mm[cc.tnum].pixsize,cc.pixsize);
	printf("\nGive latitude and longitude of image center point.\n");
	scanf("%lf%lf",&cc.clat,&cc.clon);


	if(cc.type=='a')
	{
		printf("Give standard parallels.\n");
		scanf("%lf%lf",&cc.lat1,&cc.lat2);
		if(cc.clat>0.0)
		{
			printf("Give longitude pointing down.\n");
		}
		else
		{
			printf("Give longitude pointing up.\n");
		}
		scanf("%lf",&cc.lon0);
	}
	if(cc.type=='s')	/* polar stereographic */
	{
		MU_MaxLon=180.0;
		MU_MinLon=-180.0;
		cc.scale=2.0*cc.scale*K0/cc.pixsize;
		if(cc.clat>0.0)
		{
			cc.lat1=90.0;
			printf("Give longitude pointing down.\n");
		}
		else
		{
			cc.lat1=-90.0;
			printf("Give longitude pointing up.\n");
		}
		scanf("%lf",&cc.lon0);
		cc.x0=xc;
		cc.y0=yc;
		ll_xy(&x,&y,cc.clat,cc.clon);
		cc.x0=x-xm;
		cc.y0=y-ym;
		xy_ll((double)256.5,(double)256.5,&mlat1,&mlon);
		xy_ll((double)256.5,(double)255.5,&mlat2,&mlon);
		printf("center pixel is %8.3lf km high\n",
			(mlat2-mlat1)*R*2.0*PI/360.0);
		PAUSE;
	}
	if(cc.type=='p')	/* plate carree */
	{
		MU_MaxLon=cc.clon+180.0;
		MU_MinLon=cc.clon-180.0;
		cc.scale=R*PI/180.0/cc.pixsize;	/* pixels/deg of lat */
		cc.scale2=cc.scale*cos(cc.clat*DEG_RAD); /* pixels/deg of lon */
		cc.lon0=cc.clon;
		cc.lat1=cc.clat;
		cc.x0=-xm;
		cc.y0=-ym;
	}
	if(cc.type=='m')	/* Mercator */
	{
		MU_MaxLon=cc.clon+180.0;
		MU_MinLon=cc.clon-180.0;
		cc.scale=R*PI/180.0/cc.pixsize;	/* pixels/deg of lat */
		cc.lon0=0;
		cc.lat1=0;
		cc.x0=0;
		cc.y0=0;
		ll_xy(&x,&y,cc.clat,cc.clon);
		cc.x0=x-xm;
		cc.y0=y-ym;
	}	
	if(cc.type=='l')
	{
		MU_MaxLon=cc.clon+180.0;
		MU_MinLon=cc.clon-180.0;
		if(cc.clat<0.0)sign=-1.0;
		cc.scale=R/cc.pixsize;
		cc.lon0=cc.clon;
		cc.lat1=cc.clat*sign;
		cc.n=sin(cc.lat1*DEG_RAD);
		tval=pow(tan(PI4+cc.lat1/2.0*DEG_RAD),cc.n);
		cc.F=cos(cc.lat1*DEG_RAD)*tval/cc.n;
		cc.rho0=R*cc.F/tval;
		cc.x0=xc;
		cc.y0=yc;
		ll_xy(&x,&y,cc.clat,cc.clon);
		cc.x0=x-xm;
		cc.y0=y-ym;
	}
	if(cc.type=='a')
	{
		MU_MaxLon=cc.clon+180.0;
		MU_MinLon=cc.clon-180.0;
		cc.scale=R/cc.pixsize;
		if(cc.clat<0.0)sign=-1.0;
		cc.lat1=cc.lat1*sign;
		cc.lat2=cc.lat2*sign;
		cc.n=(sin(cc.lat1*DEG_RAD)+sin(cc.lat2*DEG_RAD))/2.0;
printf("cc.n = %lf\n",cc.n);
		cc.C=cos(cc.lat1*DEG_RAD)*cos(cc.lat1*DEG_RAD)+
				2.0*cc.n*sin(cc.lat1*DEG_RAD);
printf("cc.C = %lf\n",cc.C);
		cc.rho0=R*sqrt(cc.C-2.0*cc.n*sin(cc.clat*DEG_RAD))/cc.n;
printf("cc.rho0 = %lf\n",cc.rho0);
		cc.x0=xc;
		cc.y0=yc;
		ll_xy(&x,&y,cc.clat,cc.clon);
		cc.x0=x-xm;
		cc.y0=y-ym;

	}
	if(cc.type=='S')	
	{
		MU_MaxLon=cc.clon+180;
		MU_MinLon=cc.clon-180;
		cc.scale=R*PI/180.0/cc.pixsize;	/* pixels/deg of lat */
		printf("Give longitude pointing down.\n");
		scanf("%lf",&cc.lon0);
		cc.x0=0.0;
		cc.y0=0.0;
		ll_xy(&x,&y,cc.clat,cc.clon);
		cc.x0=x-xm;
		cc.y0=y-xm;
		xy_ll((double)256.5,(double)256.5,&mlat1,&mlon);
		xy_ll((double)256.5,(double)255.5,&mlat2,&mlon);
		printf("center pixel is %8.3lf km high\n",
			(mlat2-mlat1)*R*2.0*PI/360.0);
		PAUSE;
	}

}


/********************************************************************
**
**	draw a line on the image checking for screen
**
********************************************************************* */

void cplotln(hue,x1,y1,x2,y2,val)

int hue;
int x1,y1;	/* coordinates of 1st point */
int x2,y2;	/* coordinates of 2nd point */
int val;	/* brightness of line */

{
	int   i, j, k, numx, numy;
	float x, y, dy, dx;

	if(val<0 || val >255) return;
	numx = ( x2 > x1 ) ? x2-x1 : x1-x2;
	numy = ( y2 > y1 ) ? y2-y1 : y1-y2;
	if(numx==0 && numy==0) return;
	if( numx>numy ) 
	{
		dy = (float)(y2-y1)/(float)(x2-x1);
		dx = (x2>x1) ? 1.0 : -1.0;
		i = x = x1;
		j = y = y1;
		for(k=0;k<=(numx);k++)
		{
			if(i_on_image(i,j))
				plotpt(hue,i,j,val);
			x += dx;
			y += dy*dx;
			i = x+0.5;
			j = y+0.5;
		}
	}
	else
	{
		dx = (float)(x2-x1)/(float)(y2-y1);
		dy = (y2>y1) ? 1.0 : -1.0;
		i = x = x1;
		j = y = y1;
		for(k=0;k<=(numy);k++)
		{
			if(i_on_image(i,j))
				plotpt(hue,i,j,val);
			x += dx*dy;
			y += dy;
			i = x+0.5;
			j = y+0.5;
		}
	}
}


/********************************************************************
**
**	plots a string(val,x,y,string,size) horizontal on overlay
**
********************************************************************* */

void plot_font_ovl(val,x,y,string,size,font)

int val;	         /* brightness */
int x,y;				/* location of upper right corner of 1st letter */
int size;				/* size 1 = full size  2 = 1/2 size  3 = 1/3 size */
char string[80];		/* text to plot */
int font[128][25];	/* binary font maps */

{
	int i,j,k,l,m,n;
	int mask;

	n = 0;
	while(string[n]!='\0')
	{
		for(l=0;l<25;l++)
		{
			k = string[n];
			mask = 512;
			for(m=0;m<10;m++) 
			{
				if(font[k][l]&mask) 
				{
					if(l%size == 0) 
					{
						plotpt(3,x+m/size,y+l/size,val);
					}
				}
				mask /= 2;
			}
		}
		n++;
		x += 15/size;
	}
}