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

# define MAX 100
int prec(double, double *, double *, double, double *, double *);
int inp_ra0(double *, double *, double *);
int inp_line(char [10] ,double *,double *,double *,double *, FILE *);
double dist(double,double,double,double);

main(void)
{
	char star[10];
	double vmag,bv,ra,dec;
	double ra0d,dec0dd,radius;
    double al50,dl50,al2000,dl2000;
	int al2000h,al2000m,al2000s,dl2000d,dl2000m,dl2000s,dl2000sgn;
	double pi=3.14159265358979;
	FILE *infp,*outfp;
	char name[]="E:/CATALOGS/2143A/CATALOG.DAT";
	char outname[]="D:/gspcfile.txt";

	
	/* RA0, dec0 ,radius の入力 */
	inp_ra0( &ra0d,  &dec0dd,&radius);

	/* RA0, dec0を1950年分点に変換 */
	prec(2000.0,&ra0d,&dec0dd,1950.0,&al50,&dl50);
    /* printf("1950 RA=%f  dec=%f \n",al50,dl50); */


    infp=fopen(name,"r");
	outfp=fopen(outname,"w");
	while (inp_line(star,&ra,&dec,&vmag,&bv,infp) != -1)
	/* １行入力する */
	{
		if((dist(ra,dec,al50,dl50)<radius)&&(bv>=0.5)&&(bv<=0.9)){
		/*２点間の角距離がradius以下で、0.5<=b-v等級<=0.9 のとき以下を実行 */
		/* 選んだ星の位置を2000年分点に変更 */
		prec(1950.0,&ra,&dec,2000.0,&al2000,&dl2000);
		/* h m s °’”への変換 */
		al2000=al2000/15.0;
		al2000h=(int)al2000;
		al2000=(al2000-(double)al2000h)*60.0;
		al2000m=(int)al2000;
		al2000s=(int)((al2000-(double)al2000m)*60.0);

		if(dl2000<0.0) dl2000sgn=-1;
		else dl2000sgn=1;
		dl2000=fabs(dl2000);
		dl2000d=(int)dl2000;
		dl2000=(dl2000-(double)dl2000d)*60.0;
		dl2000m=(int)dl2000;
		dl2000s=(int)((dl2000-(double)dl2000m)*60.0);
		dl2000d=dl2000d*dl2000sgn;

		printf("%s  ",star);
		printf("%5.2f  ",vmag);
		printf("%4.2f  ",bv);
		printf("%02d %02d %02d  " ,al2000h,al2000m,al2000s);
		printf("%+03d %02d %02d\n",dl2000d,dl2000m,dl2000s);

		fprintf(outfp,"%s  ",star);
		fprintf(outfp,"%5.2f  ",vmag);
		fprintf(outfp,"%4.2f  ",bv);
		fprintf(outfp,"%02d %02d %02d  " ,al2000h,al2000m,al2000s);
		fprintf(outfp,"%+03d %02d %02d\n",dl2000d,dl2000m,dl2000s);

		}
	}
    fclose(infp);
	fclose(outfp);
	return(0);
}

int prec(double t0, double *al0, double *dl0, double t, double *al, double *dl)
{
	/* 分点変換 */
	double lt0,tm,m,n,pi,al2,dl2;
	pi=3.14159265358979;

	lt0=((t+t0)/2.0-1900.0)/100.0;
	tm=(t-t0)/2.0;
	m=0.0128013+0.00000775*lt0;
	n=0.0055685-0.00000237*lt0;
	al2=*al0+(m+n*sin(*al0*pi/180.0)*tan(*dl0*pi/180.0))*tm;
	dl2=*dl0+n*cos(*al0*pi/180.0)*tm;
	*al=*al0+(m+n*sin(al2*pi/180.0)*tan(dl2*pi/180.0))*(t-t0);
	*dl=*dl0+n*cos(al2*pi/180.0)*(t-t0);
	return(0);
}

int inp_ra0(double *ra0d, double *dec0dd,double *pradius)
{
/* R.A.0 decl.0 の入力 */

	double ra0,ra0m,dec0s;
	int ra0h,dec0,dec0d,dec0m;

	printf("\nR.A.(hhmm.m)=");
	scanf("%lf",&ra0);
	/* printf("%f",ra0); */
	printf("\ndecl.(±ddmm)=");
	scanf("%d",&dec0);
	/* printf("%d",dec0); */
	printf("\nradius(dd.d)=");
	scanf("%lf",pradius);
	/* printf("%f",*pradius); */
	ra0h=(int)ra0/100 ;
	ra0m=ra0-(double)(ra0h*100) ;
	*ra0d=(double)(15*ra0h)+0.25*ra0m ;
	/* printf("\n%d h %f m     %f \n",ra0h,ra0m,*ra0d); */
	if (dec0<0) dec0s=-1.0;
		else dec0s=1.0;
	dec0=abs(dec0);
	dec0d=dec0/100;
	dec0m=dec0-dec0d*100;
	*dec0dd=(double)dec0d+(double)dec0m/60.0;
	*dec0dd=*dec0dd*dec0s ;
    /* printf("\n%d d %d m    %f \n",(int)dec0s*dec0d,dec0m,*dec0dd); */
    /*以上でra0 decl0 を 度単位に変換した */
	return(0);
}

int inp_line(char star[10],double *pra,double *pdec,
			 double *pvmag,double *pbv, FILE *infp)
{
	char buff[100],str[10];
	double ras,ded;
	int rah,ram,dem,des,i;


		if (fgets(buff,MAX,infp) != NULL)
		{
			i=1;
			sscanf(buff,"%s%lf%lf%d%d%lf%s%d%d",
			star,pvmag,pbv,&rah,&ram,&ras,str,&dem,&des);
			*pra=(double)rah+ram/60.0+ras/3600.0;
			*pra=*pra*15.0;
			ded=abs(atoi(str));
			*pdec=(double)ded+dem/60.0+des/3600.0;
			if(str[0]=='-') *pdec=-*pdec;
		}
		else
		{
			i=-1;
		}
		/* eof のときは -1 を返す */
	return(i);
}

double dist(double al1,double dl1,double al2,double dl2)
{
	double pi=3.14159265358979;
	double cd1,cd2,sd1,sd2,al21;
	double cal21,sal21,sdst,sdct,cd,sd,td,d;
	al21=al2-al1;
	cd1=cos(dl1*pi/180.0);
	cd2=cos(dl2*pi/180.0);
	sd1=sin(dl1*pi/180.0);
	sd2=sin(dl2*pi/180.0);
	cal21=cos(al21*pi/180.0);
	sal21=sin(al21*pi/180.0);
	sdst=cd2*sal21;
	sdct=cd1*sd2-sd1*cd2*cal21;
	cd  =sd1*sd2+cd1*cd2*cal21;
	sd  =sqrt(sdst*sdst+sdct*sdct);
	td  =sd/cd;
	d   =atan(td);
	if(cd<0.0) d=d+pi;
	d=d*180.0/pi;
	return(d);

}