/*
	latlon.c

from pjtest.c

	ll1.bat = cl /AL latlon.c project

*/
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "proj1.h"
#include "nomouse.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 Attribute=0;
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 clr(long lval);
long atr(long lval);

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,numfile;
	double clat,clon,cx,cy;
	FILE *fp,*fpin,*fptxt;
	char file[100],string[100];

	if(argc>1)
		strcpy(file,argv[1]);
	else
	{
		if((numfile=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,string);
	fprintf(fptxt,"%s",string);
	if(pnum>=0)
	{
		printf("Constant.proj = %c\n",Constant.proj);
		printf("Give parameters.\n\n");
		get_constants(&Constant,String,&n);
		for(k=0;k<n;k++)
			fprintf(fptxt,"%s",String[k]);
		translate(&Constant,Lat,Lon,X,Y,Tx,Ty,MAX_PTS,fpin,file,fptxt);
		fclose(fptxt);
	}
}
/**********************************************************************
**
**
**
**********************************************************************/

long clr(long lval)

{
	long color=lval%1000,sign=1;

	if(lval<0)
		color*=-1;
	return(color);
}

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

long atr(long lval)

{
	long attrib=lval/1000;

	return(attrib);
}



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

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=0;
	FILE *fpout;
	char string[100],cr=0xd;
	int val,num_pts=-1,ret,rot;
	long lval;
	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;

	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;
	}
	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("\nInput File ('%s.xyc') is of type from Digitizing Palette:\n",
			file);
		printf("     ( Output from program 'digit.c' or hand built )\n");
		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("\nInput File ('%s.xyc') is of type from Scanned Image:\n",file);
		printf("     ( Output from program 'imvis.c' )\n");
		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,
				"\nInput File ('%s.xyc') is of type from Digitizing Palette:\n",
					file);
		if(fptxt)
			fprintf(fptxt,"     ( Output from program 'digit.c' or hand built )\n");
		if(fptxt)
			fprintf(fptxt,"     ( No pixels/inch found after '*' file )\n");
 		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,
				"\nInput File ('%s.xyc') is of type from Scanned Image:\n",file);
		if(fptxt)
			fprintf(fptxt,"     ( Output from program 'imvis.c' )\n");
		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);
	printf("\nDo you want to create '%s'?  (y or n)\n",string);
	if(getch()!='y')
	{
		fclose(fpin);
		return(-1);
	}
	printf(
	"\nDo you wish to use attributes rather than colors to identify lines?");
	printf("\n  (y or n -- default is 'n' )\n");
	printf(
	"\n( Usually this is no unless the lines are elevation contours. )\n\n");
	if(getch()=='y')
		Attribute=1;
	fpout=fopen(string,"wt");
	if(!fpout)
	{
		printf("Could not open '%s' to write.\n\n",string);
		exit(0);
	}
	i=0;
	printf("----------------\n");
	if(Attribute==0)
		printf("  Num    Color\n");
	else
		printf("  Num   Attibute\n");
	printf("------- --------\n");
	while(fscanf(fpin,"%ld",&lval)==1)
	{
		if(Attribute==0)
		{
			fprintf(fpout,"%4ld",clr(lval));
			printf("%c%5d    %6ld",cr,++n,clr(lval));
		}
		else
		{
			fprintf(fpout,"%4ld",atr(lval));
			printf("%c%5d    %6ld",cr,++n,atr(lval));
		}
		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);
	fclose(fpout);
}


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

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);
}