#include <math.h>
#include <stdio.h>
main()		/* Comets program */
{
	FILE *output_file;
	double jul_date,transjul_date,q,omega,i,asc_node;
	double epochjul_date,a,e,n,w,M=0.0,semdia=0.0,mag=0.0;
	char name[16],outfile[15],c='N';
	int number=1;
	struct date{
		int year;
		int month;
		double day;
	}trans,epoch;
	double cal_jd(int y,int m,double d);

	printf("**** Comets Program ****\n");
	
		printf("Enter the name of the comet: ");
		scanf("%s",name);
		printf("Enter the date of perihelion transit (yy.mm.dd.dddd):");
		scanf("%d.%d.%lf",&trans.year,&trans.month,&trans.day);
		printf("Enter the distance at perihelion q: ");
		scanf("%lf",&q);
		printf("Enter the eccentricity e: ");
		scanf("%lf",&e);
		printf("Enter the argument of perihelion omega (1950.0): ");
		scanf("%lf",&w);
		printf("Enter the longitude of ascending node Omega (1950.0): ");
		scanf("%lf",&omega);
		printf("Enter the inclination i (1950.0): ");
		scanf("%lf",&i);
		printf("Enter the daily motion n or 0.0 if unknown: ");
		scanf("%lf",&n);
		printf("Enter the epoch (yy.mm.dd): ");
		scanf("%d.%d.%lf",&epoch.year,&epoch.month,&epoch.day);

   	        (double) transjul_date = (double) cal_jd(trans.year,
			trans.month,trans.day);
	    	(double) epochjul_date = (double) cal_jd(epoch.year,
			epoch.month,epoch.day);
			a = q/(1-e);  /* convert perihelion to mean dist */

		printf("\nComet Name: %s\n",name);
		printf("Object Number: %d\n",number);
		printf("Date of Perihelion Transit (JD): %f\n",transjul_date);	
		printf("Mean Distance: %f\n",a); /* mean distance */
		printf("Eccentricity: %f\n",e);
		printf("Argument of Latitude: %f\n",w);
		printf("Longitude of Ascending Node: %f\n",omega);
		printf("Inclindation: %f\n",i);
		printf("Epoch (JD): %f\n",epochjul_date);
	

	printf("Enter name of output file: ");
	scanf("%s",outfile);
	output_file=fopen(outfile,"w");
	
	fprintf(output_file,"%9.2f %f %f %f %f %f\n",
		transjul_date,i,omega,w,a,n);
	fprintf(output_file,"%f %f %9.2f %4.2f %4.2f %s %d\n",e,M,
		epochjul_date,mag,semdia,name,number); 

}

double cal_jd(int y,int m,double d)   /* Conversion routine */
{
	double B,jul_date;

	if (m <= 2){
		m = m + 12;
		y = y - 1;
	}
	
	y = y + 1900.0;

	B = -15.0;  /* approximation only good between 1 March 1900 and */
		 /* 28 February 2100                                  */
	
        jul_date = ((floor(365.25 * y))+(floor(30.6001*(m+1)))+(B)+
		(1720996.5)+(d));
        return((double)jul_date);    
}