/*
	project.c  2-20-90

	This is a last (I hope) attempt at standardizing the lat/lon to x/y
	and back again functions for generic map use.

	pj1.bat = cl /AL project.c vimage dispio

*/

#include <stdio.h>
#include <math.h>
#include "font.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 NUM_PROJ 16

float MaxLat=90.0,MinLat=-90.0,MaxLon=180.0,MinLon=-180.0;
unsigned char Inbuff[BUFSIZ],Inbuff2[BUFSIZ],Outbuff[BUFSIZ];

int Grey=0;
unsigned char P_Buff[12][640];

struct   Color 
{ 
	unsigned char  r, g, b; 
}lut[256],lutm[256],lutc[256],luts[256]; 
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 */
};

double S[2][4]=
{
	{1.0,-1.0,1.0,-1.0},
	{1.0,1.0,-1.0,-1.0}
};

struct projection
{
	char name[80];
	char id;
}Proj[NUM_PROJ]=
{
	{                    "Mercator",'m'},
	{         "Transverse Mercator",'u'},
	{            "Oblique Mercator",'o'},
	{       "Modified Plate Carree",'p'},
	{     "Albers Equal-Area Conic",'a'},
  	{     "Lambert Conformal Conic",'l'},
	{                       "Bonne",'b'},
	{                "Orthographic",'t'},
	{               "Stereographic",'s'},
	{                    "Gnomonic",'g'},
	{        "Vertical Perspective",'V'},
	{          "Tilted Perspective",'P'},
	{"Lambert Azimuthal Equal-Area",'L'},
	{       "Azimuthal Equidistant",'z'},
	{                  "Sinusoidal",'S'},
	{                   "Polyconic",'y'}
};

void set_lut();
int choose_projection();
int get_constants();
int ll2xy();
int xy2ll();
int graticule();
int coast_lines();
int horizon();
void compress();
int expand();



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

int auto_set_constants(cc)

struct control *cc;

{
	int i,j,k;
	char mode='p';
	char ans;
	char where[100];
	double scale,sign=1.0,tval;

	cc->scale=R/cc->scale;
	if(cc->scale<=0.0)
		return(-1);
	cc->horiz=-1.0;

	if(cc->proj=='m')
	{
		return(0);
	}
	if(cc->proj=='u')
	{
		cc->lat_0*=DEG_RAD;
		cc->lon_0*=DEG_RAD;
		cc->k0=0.9996;
		return(1);
	}
	if(cc->proj=='o')
	{
		cc->lon_p=atan2(
		(cos(cc->lat_1*DEG_RAD)*sin(cc->lat_2*DEG_RAD)*cos(cc->lon_1*DEG_RAD)-
		 sin(cc->lat_1*DEG_RAD)*cos(cc->lat_2*DEG_RAD)*cos(cc->lon_2*DEG_RAD)),
		(sin(cc->lat_1*DEG_RAD)*cos(cc->lat_2*DEG_RAD)*sin(cc->lon_2*DEG_RAD)-
		 cos(cc->lat_1*DEG_RAD)*sin(cc->lat_2*DEG_RAD)*sin(cc->lon_1*DEG_RAD)));
		cc->lat_p=atan(-cos(cc->lon_p-cc->lon_1*DEG_RAD)/tan(cc->lat_1));
		cc->lon_0=cc->lon_p+PI/2.0;
		cc->k0=0.9996;
		return(2);
	}
	if(cc->proj=='p')
	{
		if(cc->lat_0==90.0||cc->lat_0==-90.0)
			return(-1);
		cc->k1=cc->scale*PI/180.0;  /* km/deg lat */
		cc->k2=cc->scale*PI/180.0*cos(cc->lat_0*DEG_RAD);  /* km/deg lon */
	}
	if(cc->proj=='a')
	{
		if(cc->lat_0<0.0)sign=-1.0;
		cc->lat_1=cc->lat_1*sign;
		cc->lat_2=cc->lat_2*sign;
		cc->n=(sin(cc->lat_1*DEG_RAD)+sin(cc->lat_2*DEG_RAD))/2.0;
		cc->C=cos(cc->lat_1*DEG_RAD)*cos(cc->lat_1*DEG_RAD)+
				2.0*cc->n*sin(cc->lat_1*DEG_RAD);
		cc->rho_0=cc->scale*sqrt(cc->C-2.0*cc->n*sin(cc->lat_0*DEG_RAD))/cc->n;
	}
	if(cc->proj=='l')
	{
		if(cc->lat_0==0.0)return(-1);
		cc->lat_0*=DEG_RAD;
		cc->lon_0*=DEG_RAD;
		if(cc->lat_0<0.0)sign=-1.0;
		cc->lat_1=cc->lat_1*sign*DEG_RAD;
		cc->lat_2=cc->lat_2*sign*DEG_RAD;
		if(cc->lat_0==cc->lat_1&&cc->lat_0==cc->lat_2)
			cc->n=sin(cc->lat_1);
		else
			cc->n=log(cos(cc->lat_1)/cos(cc->lat_2))/
					log(tan(PI4+cc->lat_2/2.0)/tan(PI4+cc->lat_1/2.0));
		cc->F=cos(cc->lat_1)*pow(tan(PI4+cc->lat_1/2.0),cc->n)/cc->n;
		cc->rho_0=R*cc->F/pow(tan(PI4+cc->lat_0/2.0),cc->n);
	}
	if(cc->proj=='b')
	{
		cc->lon_0*=DEG_RAD;
		cc->lat_1*=DEG_RAD;
	}
	if(cc->proj=='t')
	{
		cc->lon_0*=DEG_RAD;
		cc->lat_1*=DEG_RAD;
	}
	if(cc->proj=='s')
	{
		cc->scale*=cc->k0;
		cc->lon_0*=DEG_RAD;
		cc->lat_1*=DEG_RAD;
	}
	if(cc->proj=='g')
	{
		cc->lon_0*=DEG_RAD;
		cc->lat_1*=DEG_RAD;
	}
	if(cc->proj=='V')
	{
		cc->P=cc->H/cc->scale+1.0;
		cc->horiz=cc->scale*sqrt((cc->P-1.0)/(cc->P+1.0));
		cc->lon_0*=DEG_RAD;
		cc->lat_1*=DEG_RAD;
	}
	if(cc->proj=='P')
	{
		cc->P=cc->H/cc->scale+1.0;
		cc->horiz=cc->scale*sqrt((cc->P-1.0)/(cc->P+1.0));
		cc->lon_0*=DEG_RAD;
		cc->lat_1*=DEG_RAD;
		cc->gamma*=DEG_RAD;
		cc->omega*=DEG_RAD;
		cc->H=cc->scale/(cc->P-1.0);
	}
	if(cc->proj=='L')
	{
		cc->lon_0*=DEG_RAD;
		cc->lat_1*=DEG_RAD;
	}
	if(cc->proj=='z')
	{
		cc->lon_0*=DEG_RAD;
		cc->lat_1*=DEG_RAD;
	}
	if(cc->proj=='S')
	{
		cc->lon_0*=DEG_RAD;
	}
	if(cc->proj=='y')
	{
		cc->lon_0*=DEG_RAD;
		cc->lat_0*=DEG_RAD;
	}
}



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

int expand(filename,buffer)

char *filename;
unsigned char buffer[640];

{
	int i,j,k,x,y;
	FILE *fpin;
	int xyn[3];

	fpin=fopen(filename,"rb");
	if(!fpin)return(-1);

	paint_box(3,0,P_Buff,0,0,640,480);

	while( fread((char *)xyn,sizeof(int),3,fpin)==3)
	{
		fread((char *)buffer,sizeof(char),xyn[2],fpin);
		plotrow(3,xyn[0],xyn[0]+xyn[2]-1,xyn[1],buffer);
	}
	fclose(fpin);
}




/********************************************************************
**
**	compresses an image (takes out missing data
**
********************************************************************* */

void compress(filename,buffer)

char *filename;
unsigned char buffer[640];

{
	int i,j,k,x,y;
	int *npt,xyn[3];
	int numplanes=0;
	unsigned int numpts;
	FILE *fpout;
	int x1,npr;
	long total=0;
	float perc;

	fpout = fopen(filename,"w+b");
	setbuf(fpout,Outbuff);

	if(!fpout)
	{
		printf("\n\nFailure in opening %s!!!!\n\n\n",filename);
		return;
	}

/* find first point */

	x=640;
	y=0;
	while(y<480&&x==640)
	{
		getrow(3,0,639,y,buffer);
		x=0;
		while(x<640&&buffer[x]==0)x+=1;
		y+=1;
	}
	y-=1;
	if(x>639&&y>439)
	{
		fclose(fpout);
		return;
	}

/*	process the image */

	while(y<480)
	{
		getrow(3,0,639,y,buffer);
		x=0;
		while(x<640)
		{
			x1=-1;
			while(x<640&&buffer[x]==0)x+=1;
			if(x<640)
			{
				x1=x;
				numpts=0;
				while(x<640&&buffer[x]!=0)
				{
					numpts+=1;
					x+=1;
				}
			}
			if(x1>=0)
			{
				total+=numpts;
				xyn[0]=x1;xyn[1]=y;xyn[2]=numpts;
				if((npr=fwrite((char *)xyn,sizeof(int),3,fpout))<3)
				{
					printf("\n\nOutput Disk full!!!\n\n\n");
					printf("Wrote only %d of '3'\n",npr);
					printf("HIT ANY KEY TO EXIT PROGRAM\n\n");
					getch();
					return;
				}
				if((npr=fwrite((char *)buffer+x1,sizeof(char),numpts,fpout))
					<numpts)
				{
					printf("\n\nOutput Disk full!!!\n\n\n");
					printf("Wrote only %d of %d\n",npr,numpts);
					printf("HIT ANY KEY TO EXIT PROGRAM\n\n");
					getch();
					return;
				}
			}
		}
		y+=1;
	}
	fclose(fpout);
	return;
}


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

int convert_bin(val,region,xc,yc,xs,ys,clat,clon,cc)

char *region;
int val;
int xc,yc,xs,ys;
double clat,clon;
struct control cc;

{
	int i,j,k;
	double x,y;
	int ix,iy,ixo=-1,iyo=-1;
	double dx,dy,dist;
	int secnum,numpt,numsec;
	FILE *fp,*fphdr;
	float latlon[2],mll[4],mnln,mxln,mnlt,mxlt;
	char string[40];
	int doit;
	int first_seg=-1,last_seg=-1;
	long offset;
	double cx,cy;
	int ret,reto,px,py,pxo,pyo;

	strcpy(string,region);
	strcat(string,".bin");
	fp=fopen(string,"rb");
	if(!fp)return;
	setbuf(fp,Inbuff);
	strcpy(string,region);
	strcat(string,".hdr");
	fphdr=fopen(string,"rb");
		if(!fphdr)return;
	setbuf(fphdr,Inbuff2);

	ret=ll2xy(clat,clon,&cx,&cy,cc);
	cx-=xs/2;
	cy+=ys/2;
	if(ret<0)
		return(ret);

	fread((char *)&numsec,sizeof(int),1,fphdr);
	fread((char *)&numsec,sizeof(int),1,fp);

	for(j=0;j<numsec;j++)
	{
		fread((char *)&offset,sizeof(long),1,fphdr);
		fread((char *)mll,sizeof(float),4,fphdr);
		doit=1;
		mnln=mll[0];mxln=mll[1];mnlt=mll[2];mxlt=mll[3];
		if(mnln>MaxLon||mnlt>MaxLat||mxln<MinLon||mxlt<MinLat)
		{
			doit=0;
		}

		if(doit==1)
		{
			fseek(fp,offset,SEEK_SET);
			fread((char *)&numpt,sizeof(int),1,fp);
			fread((char *)mll,sizeof(float),4,fp);

			for(i=0;i<numpt;i++)
			{
				if(kbhit()!=0)
					if(getch()=='\033')
						return(-1);
				fread((char *)latlon,sizeof(float),2,fp);
				if(doit==1)
				{
					ret=ll2xy((double)latlon[0],(double)latlon[1],&x,&y,cc);
					if(ret>=0)
					{
						ix=x-cx;
						iy=cy-y;
						if(ix>=xc&&ix<xc+xs&&iy>yc&&iy<yc+ys||
								ixo>=xc&&ixo<xc+xs&&iyo>yc&&iyo<yc+ys)
						plotpt(3,ix,iy,val);
					}

					if(i>0&&ret>0&&reto>0)
					{
						plotln(3,ix,iy,ixo,iyo,val);
					}
					reto=ret;
					ixo=ix;iyo=iy;
				}
			}
		}
	}
	fclose(fp);
	fclose(fphdr);
	return(1);
}


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

int horizon(xc,yc,xs,ys,clat,clon,cc,line,debug)

int xc,yc,xs,ys,line;
double clat,clon;
struct control cc;
int debug;

{
	int i,j,k,rs;
	double lat,lon,x,y,xo,yo,H,A;
	double cx,cy,xt,yt;

	int ret,ix,iy,px,py,pxo,pyo;
	double latmax=90.0,latmin=-90.0,lonmax=180.0,lonmin=-180.0;
	double nx,ny;

	if(cc.horiz<=1.0)return(-1);
	ret=ll2xy(clat,clon,&cx,&cy,cc);
	cx-=xs/2;
	cy+=ys/2;

	if(ret<0)
		return(ret);
	lat=cc.lat_1/DEG_RAD;lon=cc.lon_0/DEG_RAD;
	ll2xy(lat,lon,&nx,&ny,cc);
	rs=cc.horiz;	
	H=cc.horiz*cc.horiz;

if(debug==1)
	printf("cc.proj = %c H=%lf rs=%d\n",cc.proj,H,rs);

	if(cc.proj=='P')
	{
		for(k=0;k<2;k++)
		{
			for(i=0;i<=rs;i++)
			{
				xt=i;
				yt=sqrt(H-xt*xt);
				A=(yt*cos(cc.gamma)+xt*sin(cc.gamma))*
					sin(cc.omega)/cc.H+cos(cc.omega);
				x=(xt*cos(cc.gamma)-yt*sin(cc.gamma))/A;
				y=(yt*cos(cc.gamma)+xt*sin(cc.gamma))*cos(cc.omega)/A;
				if(i>0)
				{
					plotln(3,(int)((nx+xo*S[0][k])-cx),(int)(cy-(ny+yo*S[1][k])),
						(int)((nx+x*S[0][k])-cx),(int)(cy-(ny+y*S[1][k])),line);
				}
				xo=x;yo=y;
			}
			xt=cc.horiz;
			yt=0.0;
			A=(yt*cos(cc.gamma)+xt*sin(cc.gamma))*
			sin(cc.omega)/cc.H+cos(cc.omega);
			x=(xt*cos(cc.gamma)-yt*sin(cc.gamma))/A;
			y=(yt*cos(cc.gamma)+xt*sin(cc.gamma))*cos(cc.omega)/A;
			plotln(3,(int)((nx+xo*S[0][k])-cx),(int)(cy-(ny+yo*S[1][k])),
				(int)((nx+x*S[0][k])-cx),(int)(cy-(ny+y*S[1][k])),line);
		}
		return(1);
	}

	if(cc.proj=='V')
	{
		for(k=0;k<4;k++)
		{
			for(i=0;i<=rs;i++)
			{
				x=i;
				y=sqrt(H-x*x);
				px=i;
				py=y;
				if(i>0)
				{
					plotln(3,(int)((nx+pxo*S[0][k])-cx),(int)(cy-(ny+pyo*S[1][k])),
						(int)((nx+px*S[0][k])-cx),(int)(cy*S[1][k]-(ny+py)),line);
				}
				pxo=px;pyo=py;
			}
			px=rs;py=0;
			plotln(3,(int)((nx+pxo*S[0][k])-cx),(int)(cy-(ny+pyo*S[1][k])),
				(int)((nx+px*S[0][k])-cx),(int)(cy-(ny+py*S[1][k])),line);
		}
		return(1);
	}
}

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

int graticule(xc,yc,xs,ys,clat,clon,cc,line,num,spacing,debug,steps,
			constant)

int xc,yc,xs,ys,line,num;
double clat,clon;
struct control cc;
float spacing,steps;
int debug;
struct control constant;

{
	int i,j,k,rs;
	double lat,lon,x,y,H;
	double del=spacing,latp=-90.0,lonp=-90.0;
	double cx,cy;
	int ret,ix,iy,px,py,pxo,pyo;
	double latmax=90.0,latmin=-90.0,lonmax=180.0,lonmin=-180.0;
	double nx,ny;

	ret=ll2xy(clat,clon,&cx,&cy,cc);
	cx-=xs/2;
	cy+=ys/2;
	if(ret<0)
		return(ret);

	MaxLat=-100;
	MinLat= 100;
	MaxLon=-190;
	MinLon= 190;

	latp=-90.0;
	do
	{
		px=py=pxo=pyo=-999;
		lonp=lonmin;j=0;
		do
		{

if(debug==1)
	printf("1--latp = %lf  lonp = %lf\n",latp,lonp);

			if(kbhit()!=0)
				if(getch()=='\033')
					return(-1);
			ret=ll2xy(latp,lonp,&x,&y,cc);
			if(ret<0)
			{
				px=py=-999;
			}
			else
			{
				px=x-cx;
				py=cy-y;
				if(latp>MaxLat)MaxLat=latp;
				if(latp<MinLat)MinLat=latp;
				if(lonp>MaxLon)MaxLon=lonp;
				if(lonp<MinLon)MinLon=lonp;
			}

if(debug==1)
	printf("1--ret=%d cx=%lf cy=%lf x=%lf y=%lf px=%d py=%d\n",
		ret,cx,cy,x,y,px,py);

			if(px!=-999&&pxo!=-999&&(px>=xc&&px<xc+xs&&py>yc&&py<yc+ys||
					pxo>=xc&&pxo<xc+xs&&pyo>yc&&pyo<yc+ys))
			{
				plotln(3,px,py,pxo,pyo,line);
			}
			pxo=px;pyo=py;
			lonp+=del/steps;j+=1;
		}while(lonp<=lonmax);
		latp+=del;
	}while(latp<=latmax);

	latp=0.0;
	do
	{
		px=py=pxo=pyo=-999;
		lonp=lonmin;j=0;
		do
		{

if(debug==1)
	printf("2--latp = %lf  lonp = %lf\n",latp,lonp);

			if(kbhit()!=0)
				if(getch()=='\033')
					return(-1);
			ret=ll2xy(latp,lonp,&x,&y,cc);
			if(ret<0)
			{
				px=py=-999;
			}
			else
			{
				px=x-cx;
				py=cy-y;
			}

if(debug==1)
	printf("2--ret=%d cx=%lf cy=%lf x=%lf y=%lf px=%d py=%d\n",
		ret,cx,cy,x,y,px,py);

			if(px!=-999&&pxo!=-999&&(px>=xc&&px<xc+xs&&py>yc&&py<yc+ys||
					pxo>=xc&&pxo<xc+xs&&pyo>yc&&pyo<yc+ys))
			{
				plotln(3,px,py,pxo,pyo,line);
			}
			pxo=px;pyo=py;
			lonp+=del/steps;j+=1;
		}while(lonp<=lonmax);
		latp-=del;
	}while(latp>=latmin);

	lonp=0.0;
	do
	{
		px=py=pxo=pyo=-999;
		latp=latmin;i=0;
		do
		{

if(debug==1)
	printf("3--latp = %lf  lonp = %lf\n",latp,lonp);

			if(kbhit()!=0)
				if(getch()=='\033')
					return(-1);
			ret=ll2xy(latp,lonp,&x,&y,cc);
			if(ret<0)
			{
				px=py=-999;
			}
			else
			{
				px=x-cx;
				py=cy-y;
			}

if(debug==1)
	printf("3--ret=%d cx=%lf cy=%lf x=%lf y=%lf px=%d py=%d\n",
		ret,cx,cy,x,y,px,py);

			if(px!=-999&&pxo!=-999&&(px>=xc&&px<xc+xs&&py>yc&&py<yc+ys||
					pxo>=xc&&pxo<xc+xs&&pyo>yc&&pyo<yc+ys))
			{
				plotln(3,px,py,pxo,pyo,line);
			}
			pxo=px;pyo=py;
			latp+=del/steps;i+=1;
		}while(latp<=latmax);
		lonp+=del;
	}while(lonp<=lonmax);

	lonp=0.0;
	do
	{
		px=py=pxo=pyo=-999;
		latp=latmin;i=0;
		do
		{

if(debug==1)
	printf("4--latp = %lf  lonp = %lf\n",latp,lonp);

			if(kbhit()!=0)
				if(getch()=='\033')
					return(-1);
			ret=ll2xy(latp,lonp,&x,&y,constant);
			if(ret<0)
			{
				px=py=-999;
			}
			else
			{
				px=x-cx;
				py=cy-y;
			}

if(debug==1)
	printf("4--ret=%d cx=%lf cy=%lf x=%lf y=%lf px=%d py=%d\n",

		ret,cx,cy,x,y,px,py);
			if(px!=-999&&pxo!=-999&&(px>=xc&&px<xc+xs&&py>yc&&py<yc+ys||
					pxo>=xc&&pxo<xc+xs&&pyo>yc&&pyo<yc+ys))
			{
				plotln(3,px,py,pxo,pyo,line);
			}
			pxo=px;pyo=py;
			latp+=del/steps;i+=1;
		}while(latp<=latmax);
		lonp-=del;
	}while(lonp>=lonmin);

	MaxLat+=spacing;
	MinLat-=spacing;
	MaxLon+=spacing;
	MinLon-=spacing;

	return(1);
}



/********************************************************************
**
**		converts x and y in kilometers to lat and lon in degrees
**
**		returns:
**
**			0 - 13 = number of successfully processed projection
**				 -1 = no such projection
**				-99 = scale is zero
**
********************************************************************* */

int xy2ll(x,y,lat,lon,cc)

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

{
	double A,B,c,D,H,M,Q;
	double sign=1.0,rho,theta,latn,latn1,diff;

	if(cc.scale==0.0)return(-99);

	if(cc.proj=='m')
	{
		*lat=90.0-2.0*atan(exp(-y/cc.scale))/DEG_RAD;
		*lon=x/cc.scale/DEG_RAD+cc.lon_0;
		return(0);
	}
	if(cc.proj=='u')
	{
		D=y/(cc.scale*cc.k0)+cc.lat_0;
		*lat=asin(sin(D)/cosh(x/(cc.scale*cc.k0)))/DEG_RAD;
		*lon=cc.lon_0+atan(sinh(x/(cc.scale*cc.k0))/cos(D))/DEG_RAD;
		return(1);
	}
	if(cc.proj=='o')
	{
		*lat=asin(sin(cc.lat_p)*tanh(y/(cc.scale*cc.k0))+
			cos(cc.lat_p)*sin(x/(cc.scale*cc.k0))/
			cosh(y/(cc.scale*cc.k0)));
		*lat/=DEG_RAD;
		*lon=cc.lon_0+atan((sin(cc.lat_p)*sin(x/(cc.scale*cc.k0))-
			cos(cc.lat_p)*sinh(y/(cc.scale*cc.k0)))/cos(x/(cc.scale*cc.k0)));
		*lon/=DEG_RAD;
		return(2);
	}
	if(cc.proj=='p')
	{
		*lat=y/cc.k1+cc.lat_0;
		*lon=x/cc.k2+cc.lon_0;
		return(3);
	}
	if(cc.proj=='a')
	{
		if(cc.lat_0<0.0)sign=-1.0;
		y*=sign;
		rho=sqrt(x*x+(cc.rho_0-y)*(cc.rho_0-y));
		theta=atan2(x,cc.rho_0-y);
		*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.lon_0;
		return(4);
	}
	if(cc.proj=='l')
	{
		if(cc.lat_0<0.0)sign=-1.0;
		y*=sign;
		rho=sqrt(x*x+(cc.rho_0-y)*(cc.rho_0-y));
		theta=atan2(x,cc.rho_0-y);
		*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+cc.lon_0;
		*lon/=DEG_RAD;
		return(5);
	}
	if(cc.proj=='b')
	{
		if(cc.lat_0<0.0)sign=-1.0;
		rho=sqrt(x*x+
			(cc.scale*tan(1.0/cc.lat_1)-y)*
			(cc.scale*tan(1.0/cc.lat_1)-y))*sign;
		*lat=tan(1.0/cc.lat_1)+cc.lat_1-rho/cc.scale;
		*lon=cc.lon_0+
			rho*atan2(x,cc.scale*tan(1.0/cc.lat_1)-y)/
			(cc.scale*cos(*lat));
		*lat/=DEG_RAD;
		*lon/=DEG_RAD;
		return(6);
	}
	if(cc.proj=='t')
	{
		rho=x*x+y*y;
		if(rho==0.0)
		{
			*lat=cc.lat_1;
			*lon=cc.lon_0;
		}
		else
		{
			rho=sqrt(rho);
			c=asin(rho/cc.scale);
			*lat=asin(cos(c)*sin(cc.lat_1)+y*sin(c)*cos(cc.lat_1)/rho);
			if(cc.lat_1==90.0)
				*lon=cc.lon_0+atan2(x,-y);
			else if(cc.lat_1==-90.0)
				*lon=cc.lon_0+atan2(x,y);
			else
				*lon=cc.lon_0+atan2(x*sin(c),
					rho*cos(cc.lat_1)*cos(c)-y*sin(c)*sin(cc.lat_1)/rho);
		}
		*lat/=DEG_RAD;
		*lon/=DEG_RAD;
		return(7);
	}
	if(cc.proj=='s')
	{
		rho=sqrt(x*x+y*y);
		if(rho==0.0)
		{
			*lat=cc.lat_1;
			*lon=cc.lon_0;
		}
		else
		{
			c=2.0*atan(rho/(2.0*cc.scale));
			*lat=asin(cos(c)*sin(cc.lat_1)+y*sin(c)*cos(cc.lat_1)/rho);
			*lon=cc.lon_0+
				atan2(x*sin(c),rho*cos(cc.lat_1)*cos(c)-y*sin(cc.lat_1)*sin(c));
			*lat/=DEG_RAD;
			*lon/=DEG_RAD;
		}
		return(8);
	}
	if(cc.proj=='g')
	{
		rho=sqrt(x*x+y*y);
		if(rho==0.0)
		{
			*lat=cc.lat_1;
			*lon=cc.lon_0;
		}
		else
		{
			c=atan(rho/cc.scale);
			*lat=asin(cos(c)*sin(cc.lat_1)+y*sin(c)*cos(cc.lat_1)/rho);
			if(cc.lat_1==90.0)
				*lon=cc.lon_0+atan2(x,-y);
			else if(cc.lat_1==-90.0)
				*lon=cc.lon_0+atan2(x,y);
			else
				*lon=cc.lon_0+
				 atan2(x*sin(c),rho*cos(cc.lat_1)*cos(c)-y*sin(cc.lat_1)*sin(c));
			*lat/=DEG_RAD;
			*lon/=DEG_RAD;
		}
		return(9);
	}
	if(cc.proj=='V')
	{
		rho=sqrt(x*x+y*y);
		if(rho==0.0)
		{
			*lat=cc.lat_1;
			*lon=cc.lon_0;
		}
		else
		{
			c=asin((cc.P-
				sqrt(1.0-rho*rho*(cc.P+1.0)/(cc.scale*cc.scale*(cc.P-1.0)))
				)/(cc.scale*(cc.P-1.0)/rho+rho/(cc.scale*(cc.P-1.0))));
			*lat=asin(cos(c)*sin(cc.lat_1)+y*sin(c)*cos(cc.lat_1)/rho);
			*lon=cc.lon_0+
				atan2(x*sin(c),
				rho*cos(cc.lat_1)*cos(c)-y*sin(cc.lat_1)*sin(c));
			*lat/=DEG_RAD;
			*lon/=DEG_RAD;
		}
		return(10);
	}
	if(cc.proj=='P')
	{
		M=cc.H*x*cos(cc.omega)/(cc.H-y*sin(cc.omega));
		Q=cc.H*y/(cc.H-y*sin(cc.omega));
		x=M*cos(cc.gamma)+Q*sin(cc.gamma);
		y=Q*cos(cc.gamma)-M*sin(cc.gamma);
		rho=sqrt(x*x+y*y);
		if(rho==0.0)
		{
			*lat=cc.lat_1;
			*lon=cc.lon_0;
		}
		else
		{
			c=asin((cc.P-
				sqrt(1.0-rho*rho*(cc.P+1.0)/(cc.scale*cc.scale*(cc.P-1.0)))
				)/(cc.scale*(cc.P-1.0)/rho+rho/(cc.scale*(cc.P-1.0))));
			*lat=asin(cos(c)*sin(cc.lat_1)+y*sin(c)*cos(cc.lat_1)/rho);
			*lon=cc.lon_0+
				atan2(x*sin(c),
				rho*cos(cc.lat_1)*cos(c)-y*sin(cc.lat_1)*sin(c));
			*lat/=DEG_RAD;
			*lon/=DEG_RAD;
		}
		return(11);
	}
	if(cc.proj=='L')
	{
		rho=sqrt(x*x+y*y);
		if(rho==0.0)
		{
			*lat=cc.lat_1;
			*lon=cc.lon_0;
		}
		else
		{
			c=2.0*asin(rho/(2.0*cc.scale));
			*lat=asin(cos(c)*sin(cc.lat_1)+(y*sin(c)*cos(cc.lat_1)/rho));
			*lon=cc.lon_0+
				atan2(x*sin(c),rho*cos(cc.lat_1)*cos(c)-y*sin(cc.lat_1)*sin(c));
		}
		*lat/=DEG_RAD;
		*lon/=DEG_RAD;
		return(12);
	}
	if(cc.proj=='z')
	{
		rho=sqrt(x*x+y*y);
		if(rho==0.0)
		{
			*lat=cc.lat_1;
			*lon=cc.lon_0;
		}
		else
		{
			c=rho/cc.scale;
			*lat=asin(cos(c)*sin(cc.lat_1)+(y*sin(c)*cos(cc.lat_1)/rho));
			*lon=cc.lon_0+atan2(x*sin(c),
				rho*cos(cc.lat_1)*cos(c)-y*sin(cc.lat_1)*sin(c));
		}
		*lat/=DEG_RAD;
		*lon/=DEG_RAD;
		return(13);
	}
	if(cc.proj=='S')
	{
		*lat=y/cc.scale;
		if(*lat==90.0||*lat==-90.0)
			*lon=0.0;
		else
			*lon=cc.lon_0+x/(cc.scale*cos(*lat));
		*lat/=DEG_RAD;
		*lon/=DEG_RAD;
		return(14);
	}
	if(cc.proj=='y')
	{
		if(y==(-cc.scale*cc.lat_0))
		{
			*lat=0.0;
			*lon=x/cc.scale+cc.lon_0;
		}
		else
		{
			A=cc.lat_0+y/cc.scale;
			B=x*x/(cc.scale*cc.scale)+A*A;
			latn=A;
			do
			{
				latn1=latn-(A*(latn*tan(latn)+1.0)-latn-
							0.5*(latn*latn+B)*tan(latn))/
							((latn-A)/tan(latn)-1.0);
				diff=(latn-latn1)*(latn-latn1);
				latn=latn1;
			}while(diff>3.046E-18);
			*lat=latn;
			*lon=asin(x*tan(latn)/cc.scale)/sin(latn)+cc.lon_0;
		}
		*lat/=DEG_RAD;
		*lon/=DEG_RAD;
		return(15);
	}
	return(-1);
}





/********************************************************************
**
**  converts lat and lon in degrees into x and y in kilometers
**
**  returns number of projection used.  If < 0 error or back side of earth
**
**		returns:
**
**			0 - 13 = number of successfully processed projection
**				 -7 = back side of earth
**				-99 = scale is zero
**
**
********************************************************************* */

int ll2xy(lat,lon,x,y,cc)

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

{
	double A,B,c,E,k,xt,yt;
	double sign=1.0,rho,theta,tval,dlon,cosc;

	if(cc.scale==0.0)return(-99);

	if(lat>90.0||lat<-90.0)
		return(-1);

	if(cc.proj=='m')
	{
		if(lat==90.0||lat==-90.0)return(-100);
		*x=cc.scale*(lon-cc.lon_0)*DEG_RAD;
		*y=cc.scale*(log(tan((45.0+lat/2)*DEG_RAD)));
		return(0);
	}
	if(cc.proj=='u')
	{
		lat*=DEG_RAD;lon*=DEG_RAD;
		dlon=lon-cc.lon_0;
		if(dlon<=-PI2||dlon>=PI2)return(-101);
		B=cos(lat)*sin((dlon));
		if(B>=1.0||B<=-1.0)return(-101);
		*x=0.5*cc.scale*cc.k0*log((1.0+B)/(1.0-B));
		*y=cc.scale*cc.k0*(atan2(tan(lat),cos(dlon))-cc.lat_0);
		return(1);
	}
	if(cc.proj=='o')
	{
		A=sin(cc.lat_p)*sin(lat*DEG_RAD)-cos(cc.lat_p)*cos(lat*DEG_RAD)*
			sin(lon*DEG_RAD-cc.lon_0);
		if(A==1.0)return(-102);
		*x=cc.scale*cc.k0*atan((tan(lat*DEG_RAD)*cos(cc.lat_p)+
			sin(cc.lat_p)*sin(lon*DEG_RAD-cc.lon_0))/cos(lon*DEG_RAD-cc.lon_0));
		*y=cc.scale/2.0*cc.k0*log((1.0+A)/(1.0-A));
		return(2);
	}
	if(cc.proj=='p')
	{
		*x=(lon-cc.lon_0)*cc.k2;
		*y=(lat-cc.lat_0)*cc.k1;
		return(3);
	}
	if(cc.proj=='a')
	{
		if(cc.n==0.0)return(-104);
		if(cc.lat_0<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.lon_0)*DEG_RAD;
		*x=rho*sin(theta);
		*y=(cc.rho_0-rho*cos(theta))*sign;
		return(4);
	}
	if(cc.proj=='l')
	{
		lat*=DEG_RAD;lon*=DEG_RAD;
		if(cc.lat_0<0.0)sign=-1.0;
		lat*=sign;
		if((PI4+lat/2.0)<=0.0)return(-105);
		tval=pow(tan(PI4+lat/2.0),cc.n);
		rho=cc.scale*cc.F/tval;
		theta=cc.n*(lon-cc.lon_0);
		*x=rho*sin(theta);
		*y=(cc.rho_0-rho*cos(theta))*sign;
		return(5);
	}
	if(cc.proj=='b')
	{
		lat*=DEG_RAD;lon*=DEG_RAD;
		if(cc.lat_1==0.0)return(-106);
		rho=cc.scale*(tan(1.0/cc.lat_1)+cc.lat_1-lat);
		if(rho==0.0)return(-106);
		E=cc.scale*(lon-cc.lon_0)*cos(lat)/rho;
		*x=rho*sin(E);
		*y=cc.scale*tan(1.0/cc.lat_1)-rho*cos(E);
		return(6);
	}
	if(cc.proj=='t')
	{
		lat*=DEG_RAD;lon*=DEG_RAD;
		*x=cc.scale*cos(lat)*sin(lon-cc.lon_0);
		*y=cc.scale*(cos(cc.lat_1)*sin(lat)-
				sin(cc.lat_1)*cos(lat)*cos(lon-cc.lon_0));
		if(
			sin(cc.lat_1)*sin(lat)+cos(cc.lat_1)*cos(lat)*cos(lon-cc.lon_0)
			<=0.0)return(-7);  /* point on back side of earth */
		return(7);
	}
	if(cc.proj=='s')
	{
		lat*=DEG_RAD;lon*=DEG_RAD;
		tval=1.0+sin(cc.lat_1)*sin(lat)+
				cos(cc.lat_1)*cos(lat)*cos(lon-cc.lon_0);
		if(tval==0.0)return(-108);
		k=2.0/tval;
		*x=cc.scale*k*cos(lat)*sin(lon-cc.lon_0);
		*y=cc.scale*k*(cos(cc.lat_1)*sin(lat)-
				sin(cc.lat_1)*cos(lat)*cos(lon-cc.lon_0));
		return(8);
	}
	if(cc.proj=='g')
	{
		if(cc.lat_1==PI2)
		{
			if(lat<=0.0)return(-109);
			if(lat=90.0)
				*x=*y=0.0;
			else
			{
				lat*=DEG_RAD;lon*=DEG_RAD;
				dlon=lon-cc.lon_0;
				*x=cc.scale*1.0/tan(lat)*sin(dlon);
				*y=-cc.scale*1.0/tan(lat)*cos(dlon);
			}
		}
		else if(cc.lat_1==-PI2)
		{
			if(lat>=0.0)return(-109);
			if(lat=-90.0)
				*x=*y=0.0;
			else
			{
				lat*=DEG_RAD;lon*=DEG_RAD;
				dlon=lon-cc.lon_0;
				*x=-cc.scale*1.0/tan(lat)*sin(dlon);
				*y=cc.scale*1.0/tan(lat)*cos(dlon);
			}
		}
		else
		{
			lat*=DEG_RAD;lon*=DEG_RAD;
			dlon=lon-cc.lon_0;
			cosc=sin(cc.lat_1)*sin(lat)+cos(cc.lat_1)*cos(lat)*cos(dlon);
			if(cosc<=0.00001)return(-109);
			k=1/cosc;
			*x=cc.scale*k*cos(lat)*sin(dlon);
			*y=cc.scale*k*(cos(cc.lat_1)*sin(lat)-
				sin(cc.lat_1)*cos(lat)*cos(dlon));
		}
		return(9);
	}
	if(cc.proj=='V')
	{
		lat*=DEG_RAD;lon*=DEG_RAD;
		dlon=lon-cc.lon_0;
		cosc=sin(cc.lat_1)*sin(lat)+cos(cc.lat_1)*cos(lat)*cos(dlon);
		if(cosc<1.0/cc.P)return(-110);
		k=(cc.P-1.0)/(cc.P-cosc);
		*x=cc.scale*k*cos(lat)*sin(dlon);
		*y=cc.scale*k*(cos(cc.lat_1)*sin(lat)-
			sin(cc.lat_1)*cos(lat)*cos(dlon));
		return(10);
	}
	if(cc.proj=='P')
	{
		lat*=DEG_RAD;lon*=DEG_RAD;
		dlon=lon-cc.lon_0;
		cosc=sin(cc.lat_1)*sin(lat)+cos(cc.lat_1)*cos(lat)*cos(dlon);
		if(cosc<1.0/cc.P)return(-110);
		k=(cc.P-1.0)/(cc.P-cosc);
		xt=cc.scale*k*cos(lat)*sin(dlon);
		yt=cc.scale*k*(cos(cc.lat_1)*sin(lat)-
			sin(cc.lat_1)*cos(lat)*cos(dlon));
		A=(yt*cos(cc.gamma)+xt*sin(cc.gamma))*sin(cc.omega)/cc.H+cos(cc.omega);
		*x=(xt*cos(cc.gamma)-yt*sin(cc.gamma))/A;
		*y=(yt*cos(cc.gamma)+xt*sin(cc.gamma))*cos(cc.omega)/A;
		return(11);
	}
	if(cc.proj=='L')
	{
		lat*=DEG_RAD;lon*=DEG_RAD;
		dlon=lon-cc.lon_0;
		k=1.0+sin(cc.lat_1)*sin(lat)+cos(cc.lat_1)*cos(lat)*cos(dlon);
		if(k<=0.0)return(-12);
		k=sqrt(2.0/k);
		*x=cc.scale*k*cos(lat)*sin(dlon);
		*y=cc.scale*k*(cos(cc.lat_1)*sin(lat)-sin(cc.lat_1)*cos(lat)*cos(dlon));
		return(12);
	}
	if(cc.proj=='z')
	{
		lat*=DEG_RAD;lon*=DEG_RAD;
		dlon=lon-cc.lon_0;
		cosc=sin(cc.lat_1)*sin(lat)+cos(cc.lat_1)*cos(lat)*cos(dlon);
		if(cosc==-1.0)return(-113);
		if(cosc==1.0)
			*x=*y=0.0;
		else
		{
			c=acos(cosc);
			k=c/sin(c);
			*x=cc.scale*k*cos(lat)*sin(dlon);
			*y=cc.scale*k*(cos(cc.lat_1)*sin(lat)-
				sin(cc.lat_1)*cos(lat)*cos(dlon));
		}
		return(13);
	}
	if(cc.proj=='S')
	{
		lat*=DEG_RAD;lon*=DEG_RAD;
		*x=cc.scale*(lon-cc.lon_0)*cos(lat);
		*y=cc.scale*lat;
		return(14);
	}
	if(cc.proj=='y')
	{
		lat*=DEG_RAD;lon*=DEG_RAD;
		if(lat==0.0)
		{
			*x=cc.scale*(lon-cc.lon_0);
			*y=-cc.scale*cc.lat_0;
		}
		else
		{
			E=(lon-cc.lon_0)*sin(lat);
			*x=cc.scale*sin(E)/tan(lat);
			*y=cc.scale*(lat-cc.lat_0+(1.0-cos(E))/tan(lat));
		}
		return(15);
	}
	return(-1);
}





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

int get_constants(cc)

struct control *cc;

{
	int i,j,k;
	char mode='p';
	char ans;
	char where[100];
	double scale,sign=1.0,tval;

	if(cc->proj=='m')strcpy(where,"at the equator");
	if(cc->proj=='u')strcpy(where,"along the central meridian");
	if(cc->proj=='o')strcpy(where,"along the great circle arc");
	if(cc->proj=='p')strcpy(where,"at the map or image center");
	if(cc->proj=='a')strcpy(where,"at the map or image center");
	if(cc->proj=='l')strcpy(where,"at the map or image center");
	if(cc->proj=='b')strcpy(where,"at the map or image center");
	if(cc->proj=='t')strcpy(where,"at the map or image center");
	if(cc->proj=='s')strcpy(where,"at the pole");
	if(cc->proj=='g')strcpy(where,"at the map or image center");
	if(cc->proj=='V')strcpy(where,"at the map or image center");
	if(cc->proj=='P')strcpy(where,"at the map or image center");
	if(cc->proj=='L')strcpy(where,"at the map or image center");
	if(cc->proj=='z')strcpy(where,"at the map or image center");
	if(cc->proj=='S')strcpy(where,"at the map or image center");
	if(cc->proj=='y')strcpy(where,"at the map or image origin");

	printf("Enter scale of map ( 1:______) or pixel size %s.\n",where);
	scanf("%lf",&scale);
	cc->scale=R/scale;
	if(cc->scale<=0.0)
		return(-1);
	cc->horiz=-1.0;

	if(cc->proj=='m')
	{
		printf("Give longitude (degrees) of central meridian of map\n");
		scanf("%lf",&cc->lon_0);
		return(0);
	}
	if(cc->proj=='u')
	{
		printf("Give central latitude (degrees) of map\n");
		scanf("%lf",&cc->lat_0);
		cc->lat_0*=DEG_RAD;
		printf("Give longitude (degrees) of central meridian of map\n");
		scanf("%lf",&cc->lon_0);
		cc->lon_0*=DEG_RAD;
		cc->k0=0.9996;
		return(1);
	}
	if(cc->proj=='o')
	{
		printf("Give latitude and longitude (degrees) of first point.\n");
		scanf("%lf%lf",&cc->lat_1,&cc->lon_1);
		printf("Give latitude and longitude (degrees) of second point.\n");
		scanf("%lf%lf",&cc->lat_2,&cc->lon_2);
		cc->lon_p=atan2(
		(cos(cc->lat_1*DEG_RAD)*sin(cc->lat_2*DEG_RAD)*cos(cc->lon_1*DEG_RAD)-
		 sin(cc->lat_1*DEG_RAD)*cos(cc->lat_2*DEG_RAD)*cos(cc->lon_2*DEG_RAD)),
		(sin(cc->lat_1*DEG_RAD)*cos(cc->lat_2*DEG_RAD)*sin(cc->lon_2*DEG_RAD)-
		 cos(cc->lat_1*DEG_RAD)*sin(cc->lat_2*DEG_RAD)*sin(cc->lon_1*DEG_RAD)));
		cc->lat_p=atan(-cos(cc->lon_p-cc->lon_1*DEG_RAD)/tan(cc->lat_1));
		cc->lon_0=cc->lon_p+PI/2.0;
		cc->k0=0.9996;
		return(2);
	}
	if(cc->proj=='p')
	{
		printf("Give latitude and longitude (degrees) of origin.\n");
		scanf("%lf%lf",&cc->lat_0,&cc->lon_0);
		if(cc->lat_0==90.0||cc->lat_0==-90.0)
			return(-1);
		cc->k1=cc->scale*PI/180.0;  /* km/deg lat */
		cc->k2=cc->scale*PI/180.0*cos(cc->lat_0*DEG_RAD);  /* km/deg lon */
	}
	if(cc->proj=='a')
	{
		printf("Give latitude and longitude (degrees) of origin.\n");
		scanf("%lf%lf",&cc->lat_0,&cc->lon_0);
		printf("Give latitude of 2 standard parallels.\n");
		scanf("%lf%lf",&cc->lat_1,&cc->lat_2);
		if(cc->lat_0<0.0)sign=-1.0;
		cc->lat_1=cc->lat_1*sign;
		cc->lat_2=cc->lat_2*sign;
		cc->n=(sin(cc->lat_1*DEG_RAD)+sin(cc->lat_2*DEG_RAD))/2.0;
		cc->C=cos(cc->lat_1*DEG_RAD)*cos(cc->lat_1*DEG_RAD)+
				2.0*cc->n*sin(cc->lat_1*DEG_RAD);
		cc->rho_0=cc->scale*sqrt(cc->C-2.0*cc->n*sin(cc->lat_0*DEG_RAD))/cc->n;
	}
	if(cc->proj=='l')
	{
		printf("Give latitude and longitude (degrees) of origin.\n");
		scanf("%lf%lf",&cc->lat_0,&cc->lon_0);
		if(cc->lat_0==0.0)return(-1);
		cc->lat_0*=DEG_RAD;
		cc->lon_0*=DEG_RAD;
		printf("Give latitude of 2 standard parallels.\n");
		scanf("%lf%lf",&cc->lat_1,&cc->lat_2);
		if(cc->lat_0<0.0)sign=-1.0;
		cc->lat_1=cc->lat_1*sign*DEG_RAD;
		cc->lat_2=cc->lat_2*sign*DEG_RAD;
		if(cc->lat_0==cc->lat_1&&cc->lat_0==cc->lat_2)
			cc->n=sin(cc->lat_1);
		else
			cc->n=log(cos(cc->lat_1)/cos(cc->lat_2))/
					log(tan(PI4+cc->lat_2/2.0)/tan(PI4+cc->lat_1/2.0));
		cc->F=cos(cc->lat_1)*pow(tan(PI4+cc->lat_1/2.0),cc->n)/cc->n;
		cc->rho_0=R*cc->F/pow(tan(PI4+cc->lat_0/2.0),cc->n);
	}
	if(cc->proj=='b')
	{
		printf("Give latitude of standard parallel.\n");
		scanf("%lf",&cc->lat_1);
		printf("Give longitude (degrees) of central meridian.\n");
		scanf("%lf",&cc->lon_0);
		cc->lon_0*=DEG_RAD;
		cc->lat_1*=DEG_RAD;
	}
	if(cc->proj=='t')
	{
		printf("Give latitude of standard parallel.\n");
		scanf("%lf",&cc->lat_1);
		printf("Give longitude (degrees) of central meridian.\n");
		scanf("%lf",&cc->lon_0);
		cc->lon_0*=DEG_RAD;
		cc->lat_1*=DEG_RAD;
	}
	if(cc->proj=='s')
	{
		printf("Give latitude of standard parallel.\n");
		scanf("%lf",&cc->lat_1);
		printf("Give longitude (degrees) of central meridian.\n");
		scanf("%lf",&cc->lon_0);
		printf("GIve scale factor (usually 1).\n");
		scanf("%lf",&cc->k0);
		cc->scale*=cc->k0;
		cc->lon_0*=DEG_RAD;
		cc->lat_1*=DEG_RAD;
	}
	if(cc->proj=='g')
	{
		printf("Give latitude of standard parallel.\n");
		scanf("%lf",&cc->lat_1);
		printf("Give longitude (degrees) of central meridian.\n");
		scanf("%lf",&cc->lon_0);
		cc->lon_0*=DEG_RAD;
		cc->lat_1*=DEG_RAD;
	}
	if(cc->proj=='V')
	{
		printf("Give latitude of standard parallel.\n");
		scanf("%lf",&cc->lat_1);
		printf("Give longitude (degrees) of central meridian.\n");
		scanf("%lf",&cc->lon_0);
		printf("GIve altitude of view in kilometers.\n");
		scanf("%lf",&cc->H);
		cc->P=cc->H/cc->scale+1.0;
		cc->horiz=cc->scale*sqrt((cc->P-1.0)/(cc->P+1.0));
		cc->lon_0*=DEG_RAD;
		cc->lat_1*=DEG_RAD;
	}
	if(cc->proj=='P')
	{
		printf("Give latitude of standard parallel.\n");
		scanf("%lf",&cc->lat_1);
		printf("Give longitude (degrees) of central meridian.\n");
		scanf("%lf",&cc->lon_0);
		printf("GIve altitude of view in kilometers.\n");
		scanf("%lf",&cc->H);
		cc->P=cc->H/cc->scale+1.0;
		cc->horiz=cc->scale*sqrt((cc->P-1.0)/(cc->P+1.0));
		cc->lon_0*=DEG_RAD;
		cc->lat_1*=DEG_RAD;
		printf("Give viewing azimuth (N=0,E=90,S=180,W=270).\n");
		scanf("%lf",&cc->gamma);
		printf("Give viewing angle up from straight down.\n");
		scanf("%lf",&cc->omega);
		cc->gamma*=DEG_RAD;
		cc->omega*=DEG_RAD;
		cc->H=cc->scale/(cc->P-1.0);
	}
	if(cc->proj=='L')
	{
		printf("Give latitude of standard parallel.\n");
		scanf("%lf",&cc->lat_1);
		printf("Give longitude (degrees) of central meridian.\n");
		scanf("%lf",&cc->lon_0);
		cc->lon_0*=DEG_RAD;
		cc->lat_1*=DEG_RAD;
	}
	if(cc->proj=='z')
	{
		printf("Give latitude of standard parallel.\n");
		scanf("%lf",&cc->lat_1);
		printf("Give longitude (degrees) of central meridian.\n");
		scanf("%lf",&cc->lon_0);
		cc->lon_0*=DEG_RAD;
		cc->lat_1*=DEG_RAD;
	}
	if(cc->proj=='S')
	{
		printf("Give longitude (degrees) of central meridian.\n");
		scanf("%lf",&cc->lon_0);
		cc->lon_0*=DEG_RAD;
	}
	if(cc->proj=='y')
	{
		printf("Give longitude (degrees) of central meridian.\n");
		scanf("%lf",&cc->lon_0);
		printf("Give latitude of origin ( usually = 0.0).\n");
		scanf("%lf",&cc->lat_0);
		cc->lon_0*=DEG_RAD;
		cc->lat_0*=DEG_RAD;
	}
}





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

int choose_projection(cc)

struct control *cc;

{
	int i,num=-1;
	char ans;

	printf("\nChoose projection type\n\n");
	for(i=0;i<NUM_PROJ;i++)
		printf("     %c -- %s\n",Proj[i].id,Proj[i].name);
	ans=getche();
	printf("\n");
	for(i=0;i<NUM_PROJ;i++)
		if(ans==Proj[i].id)
			num=i;
	cc->proj=ans;
	cc->pnum=num;

	return(num);
}



/********************************************************************
**
**	type  g -- grey all lut[] files
**			G -- grey only lut[]
**			i -- initial header map
**			r -- reset lut[] to lutc[]
**
********************************************************************* */

void set_lut(type)

char type;

{
	int k,hue,brite,j;
	int color[3];
	double dbrite,dval;
	float val,minval;

	if(type=='c')
	{
		Grey=0;
		for(k=15;k<256;k++) 
		{ 
			lutc[k].r=lut[k].r=k; 
			lutc[k].g=lut[k].g=k; 
			lutc[k].b=lut[k].b=k; 
		} 
		lut[ 0].r=  0;lut[ 0].g=  0;lut[ 0].b=  0;
		lut[ 1].r= 46;lut[ 1].g= 46;lut[ 1].b= 46;
		lut[ 2].r=128;lut[ 2].g=128;lut[ 2].b=128;
		lut[ 3].r=255;lut[ 3].g=255;lut[ 3].b=255;
		lut[ 4].r=255;lut[ 4].g=  0;lut[ 4].b=  0;
		lut[ 5].r=255;lut[ 5].g=128;lut[ 5].b=  0;
		lut[ 6].r=255;lut[ 6].g=255;lut[ 6].b=  0;
		lut[ 7].r=  0;lut[ 7].g=255;lut[ 7].b=  0;
		lut[ 8].r=  0;lut[ 8].g=255;lut[ 8].b=255;
		lut[ 9].r= 64;lut[ 9].g= 64;lut[ 9].b=255;
		lut[10].r=255;lut[10].g=  0;lut[10].b=255;
		lut[11].r=192;lut[11].g= 96;lut[11].b= 48;
		lut[12].r= 96;lut[12].g= 28;lut[12].b= 14;
		lut[13].r=128;lut[13].g=128;lut[13].b=255;
		lut[14].r=105;lut[14].g= 90;lut[14].b= 75;
		lut[15].r=  0;lut[15].g= 20;lut[15].b= 60;
	}
	if(type=='g')
	{
		for(k=0;k<256;k++) 
		{ 
			lutc[k].r=lut[k].r=k; 
			lutc[k].g=lut[k].g=k; 
			lutc[k].b=lut[k].b=k; 
		} 
		Grey=1;
	}
	if(type=='G')
	{
		for(k=0;k<256;k++) 
		{ 
			lut[k].r=k; 
			lut[k].g=k; 
			lut[k].b=k; 
		} 
		Grey=1;
	}
	if(type=='i')
	{
		for(k=0;k<256;k++) 
		{ 
			lutc[k].r=lut[k].r=k; 
			lutc[k].g=lut[k].g=k; 
			lutc[k].b=lut[k].b=k; 
		} 
		lut[ 64].r=  0;lut[ 64].g=  0;lut[ 64].b=128;
		lut[127].r=  0;lut[127].g=100;lut[127].b=  0;
		lut[191].r=255;lut[191].g=  0;lut[191].b=  0;
		lut[192].r=255;lut[192].g=128;lut[192].b=  0;
		lut[193].r=255;lut[193].g=255;lut[193].b=  0;
		lut[194].r=  0;lut[194].g=255;lut[194].b=  0;
		lut[195].r=  0;lut[195].g=255;lut[195].b=255;
		lut[196].r= 20;lut[196].g= 55;lut[196].b=255;
		lut[253].r=100;lut[253].g=  1;lut[253].b=  1;
		lut[254].r=100;lut[254].g=  1;lut[254].b=  1;
	}
	if(type=='r')
	{
		for(k=0;k<256;k++) 
		{ 
			lut[k].r=lutc[k].r; 
			lut[k].g=lutc[k].g; 
			lut[k].b=lutc[k].b;
		} 
	}
	WritePalette(lut); 
}