
use Math::Trig;


# Julian Date (JD) is the number of days from 12:00 UT on 1 Jan 4713 BC.
# conversion from textbook Murray and Dermott.
sub calendar2julian($$$$) {
	my($y0,$m0,$d0,$ut0);
	my($y1,$m1,$b);
	my($ans);

	$y0=1.0*$_[0];
	$m0=1.0*$_[1];
	$d0=1.0*$_[2];
	$ut0=1.0*$_[3];

	if ($m0 > 2) {
		$y1=$y0;
		$m1=$m0;
	} else {
		$y1=$y0-1;
		$m1=$m0+12;
	}

	if (($y0 <= 1582) && ($m0 <= 10) && ($d0 <=4)) {
		$b=-2;
	} else {
		$b=floor($y1/400)-floor($y1/100);
	}
# print  "$y1, $m1, $b\n";
	$ans=floor(365.25*$y1)+floor(30.6001*($m1+1))+$b+1720996.5+$d0+$ut0/24.0;


	return $ans;
}

# Julian Date (JD) is the number of days from 12:00 UT on 1 Jan 4713 BC.
sub julian2calendar($) {
	my($a, $b, $c, $d, $e, $f);
	my($jd);

	$jd=1.0*$_[0];

	$a=floor($jd+0.5);
	if($a < 2299161) {
		$c=$a+1524;
	} else {
		$b=floor(($a-1867216.25)/36524.25);
		$c=$a+$b-floor($b/4)+1525;
	}
	$d=floor(($c-122.1)/365.25);
	$e=floor(365.25*$d);
	$f=floor(($c-$e)/30.6001);

# print "$a, $b, $c, $d, $e, $f\n";
	$jd_day=$c-$e-floor(30.6001*$f)+frac($jd+0.5);
	$jd_month=$f-1-12*floor($f/14);
	$jd_year=$d-4715-floor((7.0+$jd_month)/10);
	
	return 13;
}

# convert the eccentric anomaly (E) to true anomaly (f).
sub e2f($$) {
	my($e,$ecc);
	$e=$_[0];
	$ecc=$_[1];
	
	return 2.0*atan(((1.0+$ecc)/(1.0-$ecc))**0.5*tan($e/2.0));
}

# Quick Newton-Raphson solution to the Kepler Equation.
sub solve_kepler($$) {
	my($m,$e,$ecc);
	my($f,$g, $i);
	my($e0, $diff);
	
	$m=$_[0];
	$ecc=$_[1];
	
	if ($m < $pi) {
		$e0=$pi/2.0;
	} else { 
		$e0=3.0*$pi/2.0;
	}
	
	do {
#print "$i	$e0	";
		$f=$e0-$ecc*sin($e0)-$m;
		$g=1.0-$ecc*cos($e0);
		$diff=$f/$g;
		$e0=$e0-$diff;
#printf "%g	$e0\n",$f/$g;
	} until (abs($diff) < 1e-14);
	
	return $e0;
}

sub xyz2sph ($$$) {
	my($x, $y, $z);
	
	$x=$_[0];
	$y=$_[1];
	$z=$_[2];
	
# print "$x, $y, $z\n";
	
	$sph_r=sqrt($x*$x + $y*$y + $z*$z);
	($sph_r==0) && return 0;
	$sph_lat=asin($z/$sph_r);
	$sph_long=atan($y/$x);
	if ($x<0) { $sph_long+=pi(); }
	if ($sph_long < 0) { $sph_long+=2.0*pi(); };
	
	return $sph_r;
}

sub sph2xyz ($$$) {
	my($lat, $lng, $r, $rcoslat);
	
	$lat=$_[0];
	$lng=$_[1];
	$r=$_[2];
	
# printf "%g, %g, $r\n",rad2deg($lat),rad2deg($lng);

	$cat_z=$r*sin($lat);
	$rcoslat=$r*cos($lat);
	$cat_x=$rcoslat*cos($lng);
	$cat_y=$rcoslat*sin($lng);
	
	return $sph_r;
}

sub select_time() {
	my($ans, $ut, $day);
	my($i, $line);
	local(*fp);
	
	$ans=read_input("Enter Julian date (jd), or calendar date (cal), or (list)?","cal");
	
	if ($ans eq "jd") {
		$time_jd=read_input("Enter Julian Date ",$time_jd);
		julian2calendar($time_jd);
		$time_year=$jd_year;
		$time_month=$jd_month;
		$time_day=$jd_day;
		print "The calendar date is $time_year, $time_month, $time_day\n";
	}elsif ($ans eq "cal") {
		$time_year=read_input("Enter year ",$time_year);
		$time_month=read_input("Enter month ",$time_month);
		$time_day=read_input("Enter day (use fraction to indicate time of day)",$time_day);
		$day=floor($time_day);
		$ut=24.0*($time_day-$day);
		$time_jd=calendar2julian($time_year,$time_month,$day,$ut);
		print "The Julian Date is $time_jd\n";
	} elsif ($ans eq "list") {
		$ans=read_input("What file has the list of JD?","$list_jd");
		if (open(fp,"$ans")) {
			print "Found file $ans, first set of lines...\n";
			$i=0;
			while ($line=<fp>) {
				print $line;
				++$i;
				($i==20) && close(fp);
			}
			$time_jd=$time_year=$time_month="list";
			$list_jd=$ans;
		} else {
			print "didn't find $ans\n";
		}
	} else { 
		print "\n*** Sorry, didn't understand what date mode was requested\n\n";
	}
	
	return 1;
}




