/*
	albers.c

	converts x,y to lat,lon
	or lat,lon to x,y


	To use a ____.bat file to run you can name the output file
	and either make a new file or append data to an existing one.

	albers s10_str  str10_13  /n
	albers s11_str  str10_13  /a

	This will create a new file 'tr10_13.raw' and put both
	's10_str.txt' and 's11_str.txt' into it.

*/
#include <stdio.h>
#include <math.h>
#include <string.h>

#define PI		3.1415927
#define PI4		0.7853981
#define DEG_RAD 0.01745329251				/* converts degrees to radians */
#define NUM_PROJ 10
#define SKIP 	while(fgetc(fpin)!=10)


struct proj
{
	char name[20];
	int year;
	double a,b;
	char remark[100];
}Ellip[NUM_PROJ]=
{
	{"GRS 80",1980,6378137.0,6356752.3,"Newly adopted"},
	{"WGS 72",1972,6378135.0,6356750.5,"NASA, DoD, oil Co."},
	{"Australian",1965,6378160.0,6356774.7,"Australia"},
	{"Krasovsky",1940,6378245.0,6356863.0,"Soviet Union"},
	{"International",1924,6378388.0,6356911.9,"Remainer of the World"},
	{"Hayford",1924,6378388.0,6356911.9,"Remainer of the World"},
	{"Clarke",1880,6378249.1,6356514.9,"Most of Africa; France"},
	{"Clarke",1866,6378206.4,6356583.8,"North America; Philippines"},
	{"Airy",1830,6377563.4,6356256.9,"Great Britain"},
	{"Everest",1830,6377276.3,6356075.4,"India;Burma;Pak.;Afgan.;Thailand;etc."}
};
int Proj;
double e,a,e2,n,C,rho0,q_max;

int set_p(void);
double Lat0,Lat1,Lat2,Lon0;
int set_const();
int axy_ll_e();
int all_xy_e();
int axy_ll_e();

/***********************************************************************
**
**
**
**
************************************************************************/

int test_albers()


{
	int i,j,k;
	double x,y,lat,lon;
	FILE *fpin,*fpout;
	char file[100],infile[100],outfile[100],string[100],ans;
	int n[5],attrib1,attrib2;

	Proj=7;
	printf("Proj = %d\n\n",Proj);

	Lat0=23.0;
	Lat1=29.5;
	Lat2=45.5;
	Lon0=-96.0;
	set_const(Proj,Lat0,Lat1,Lat2,Lon0);	
	printf("a = %lf e = %20.17lf e**2 = %20.17lf\n",a,e,e2);

	do
	{
		printf("   1 --  x/y - lat/lon\n");
		printf("   2 --  lat/lon to x/y\n");
		printf("\n   x -- exit program\n\n");
		ans=getch();
		if(ans=='1')
		{
			do
			{
				printf("Give x and y (9999 to end)\n");
				scanf("%lf",&x);
				if(x!=9999.0)
				{
					scanf("%lf",&y);
					axy_ll_e(x,y,&lat,&lon);
					printf("lat = %10.6lf  lon = %11.6lf\n",lat,lon);
				}
			}while(x!=9999.0);
		}
		if(ans=='2')
		{
			do
			{
				printf("Give lat and lon (9999 to end)\n");
				scanf("%lf",&lat);
				if(lat!=9999.0)
				{
					scanf("%lf",&lon);
					all_xy_e(&x,&y,lat,lon);
					printf("x = %10.4lf  y = %10.4lf\n",x,y);
				}
			}while(lat!=9999.0);
		}
	}while(ans!='x');

}

/********************************************************************
**
**
**
********************************************************************* */

int axy_ll_e(x,y,lat,lon)

double x,y,*lat,*lon;

{
	double q,rho,theta;
	double ry=rho0-y,err,sign=1.0,tlat,sin_tlat;
	int time=0,maxtime=100;

	if(n<0.0)sign=-1.0;
	theta=atan2(x*sign,ry*sign);
	rho=sqrt(x*x+ry*ry);
	q=(C-rho*rho*n*n/(a*a))/n;

	if(q>2.0||q<-2.0)
	{
		return(-1);
	}
	tlat=asin(q/2.0);

	if(q<-q_max||q>q_max)
	{
		return(-1);
	}
	else if(q==q_max)
		*lat=90.0*DEG_RAD;
	else if(q==-q_max)
		*lat=-90.0*DEG_RAD;
	else do
	{
		sin_tlat=sin(tlat);
		*lat=tlat+(1.0-e2*sin_tlat*sin_tlat)*(1.0-e2*sin_tlat*sin_tlat)/
			(2.0*cos(tlat))*
			(q/(1.0-e2)-sin(tlat)/(1.0-e2*sin_tlat*sin_tlat)+1/(2.0*e)*
			log((1.0-e*sin_tlat)/(1.0+e*sin_tlat)));
		err=(*lat-tlat)*(*lat-tlat);
		tlat=*lat;
	}while(err>0.000000001&&time++<maxtime);

	*lon=Lon0+theta/n;
	*lon/=DEG_RAD;
	*lat/=DEG_RAD;

	return(1.0);
}




/********************************************************************
**
**
**
********************************************************************* */

int all_xy_e(x,y,lat,lon)

double *x,*y,lat,lon;

{
	double q,m,rho,theta;

	if(lat>90.0||lat<-90.0)
		return(-1);

	lat*=DEG_RAD;
	lon*=DEG_RAD;

	q=(1.0-e*e)*(sin(lat)/(1.0-e*e*sin(lat)*sin(lat))-1.0/(2.0*e)*
		(log((1.0-e*sin(lat))/(1.0+e*sin(lat)))));
	m=cos(lat)/sqrt(1.0-e*e*sin(lat)*sin(lat));
	rho=a*sqrt(C-n*q)/n;
	theta=n*(lon-Lon0);
	*x=rho*sin(theta);
	*y=rho0-rho*cos(theta);

	return(1);
}


/********************************************************************
**
**	input projection number 
**	plus coordinates in degrees.
**
********************************************************************* */

int set_const(p,lat0,lat1,lat2,lon0)

int p;
double lat0,lat1,lat2,lon0;

{
	double m1,m2,q0,q1,q2;

	lat0*=DEG_RAD;
	lat1*=DEG_RAD;
	lat2*=DEG_RAD;
	lon0*=DEG_RAD;
	Lat0=lat0;
	Lat1=lat1;
	Lat2=lat2;
	Lon0=lon0;

	a=Ellip[p].a;
	e2=1.0-Ellip[p].b*Ellip[p].b/(Ellip[p].a*Ellip[p].a);
	e=sqrt(e2);
	q0=(1.0-e*e)*(sin(lat0)/(1.0-e*e*sin(lat0)*sin(lat0))-1.0/(2.0*e)*
		(log((1.0-e*sin(lat0))/(1.0+e*sin(lat0)))));
	q1=(1.0-e*e)*(sin(lat1)/(1.0-e*e*sin(lat1)*sin(lat1))-1.0/(2.0*e)*
		(log((1.0-e*sin(lat1))/(1.0+e*sin(lat1)))));
	q2=(1.0-e*e)*(sin(lat2)/(1.0-e*e*sin(lat2)*sin(lat2))-1.0/(2.0*e)*
		(log((1.0-e*sin(lat2))/(1.0+e*sin(lat2)))));
	m1=cos(lat1)/sqrt(1.0-e*e*sin(lat1)*sin(lat1));
	m2=cos(lat2)/sqrt(1.0-e*e*sin(lat2)*sin(lat2));
	if(lat1==lat2)
		n=sin(lat1);
	else
		n=(m1*m1-m2*m2)/(q2-q1);

	C=m1*m1+n*q1;
	rho0=a*sqrt(C-n*q0)/n;
	q_max=1.0-((1.0-e2)/(2.0*e))*log((1.0-e)/(1.0+e));
}


/********************************************************************
**
**
**
********************************************************************* */

int set_p()

{
	int i,j,k,p;

	do
	{
		for(i=0;i<NUM_PROJ;i++)
			printf("%2d%15s %4d %9.1lf %9.1lf %s\n",
				i+1,Ellip[i].name,Ellip[i].year,Ellip[i].a,Ellip[i].b,
				Ellip[i].remark);
		printf("\nEnter number of choice\n");
		scanf("%d",&p);
	}while(p<1||p>NUM_PROJ);
	return(p-1);
}







             