#!/usr/bin/perl
# by Alex Forbes and summitlake.com 2003
# noncommercial use OK so long as credits remain as supplied.
# this is the subroutine library called by the unpublished Perl program:
# "Astrophotography Exposure Tables".
# library name: astro-lib.pl
#==============================================================================
sub magnification { #p. 78 Covington
	# solve for M, projection magnification
	# s2 = distance from objective to film plane, i.e. 75mm
	# F2 = eyepiece focal length, i.e. 18mm eyepiece
	my ($s2, $F2) = @_;
	my $M = ($s2 - $F2)/$F2;
	return ( &round_to($M,3) );
}

sub magnification_afocal { #p. 76 Covington
	# solve for M, afocal magnification
	# F2 = eyepiece focal length, i.e. 18mm eyepiece
	# FC = camera lens focal length in mm (if any)
	my ($FC, $F2) = @_;
	$FC = $F2 unless ($FC);
	my $M = $FC/$F2;
	return ( &round_to($M,3) );
}

sub reciprocity {
	#correct for film reciprocity failure for t GE 1 second
	my ($t, $p) = @_;
	return ((($t+1)**(1/$p)) - 1);
}

sub system_focal { #p. 78 Covington
	# after solving for M, find overall focal length
	# F1 = telescope focal length, i.e. 1200 mm
	my ($F1, $M) = @_;
	my $F = $F1 * $M;
	return ( &round_to($F,0) );
}

sub system_fratio { #p. 78 Covington
	# after solving for M, find overall f-ratio
	# f1 = telescope f-ratio, i.e. the 8 in f/8
	my ($f1, $M) = @_;
	my $f = $f1 * $M;
	return ( &round_to($f,1) );
}

sub calc_shutter_output {
	my ($film_speed, $fratio, $B) = @_;
	print "film speed is: ", $film_speed, "<br>\n";
	print "f-ratio is: ", $fratio, "<br>\n";
	print "B is: ", $B, "<br>\n";
	print "<p>Then exposure is <br>\n";
	$t = &exposure($fratio, $film_speed, $B, 1);
	if ($t > 1) {
		#print $t, " seconds as a decimal number<br>\n";
		print (&round_to($t,0)); print " seconds as closest integer<br>\n";
	}
	else {
	$tdivisor = &divisor($t);
	print "1/",$tdivisor, " second, as a fraction<br>\n";
	print "1/",(&closest_divisor($tdivisor)), " second, as closest shutter speed<br>\n";
 	}
}

sub return_shutter_output {
	my ($film_speed, $fratio, $B, $filter_factor) = @_;
	$filter_factor = 1 unless ($filter_factor);
	$t = &exposure($fratio, $film_speed, $B, $filter_factor);
	
	if ( ($t >= 1) && $input{'is_film'} ) {
		$t = &reciprocity($t, $schwarzschild);
	}
	
	return ">> &nbsp;" if ($t > $max_secs);
	if ($t >= 1) {
		my $rt = &round_to($t,0);
		return 1 if ($rt == 1);
		return $rt;
	}
	else {
		my $td = &closest_divisor(&divisor($t));
		return 1 if ($td ==1);
		return "<< &nbsp;" if ($td == 9999);
		return "1/".$td;
	}
}

sub closest_divisor { # no good for speed > 1 second
	# for t <= 1 second, convert 1/1051 to 1/1000
	my $t = shift; my $i; my $mid; my $count;
	return $t if ($t < 1);
	my @speeds = (1,2,4,8,15,30,60,125,250,500,1000,2000,4000,8000,9999);
	$count = scalar (@speeds);
	return 9999 if ($t > $speeds[count-2]); # limit is 1/8000
	
	for ( $i = 0; $i < $count; $i++) {
		next if ($t >= $speeds[$i]);
		$mid = ($speeds[$i] + $speeds[$i-1])/2;
		#print "<br>test: mid is $mid, t is $t, i is $i, speeds is $speeds[$i]<br>\n";
		if ($t < $mid) {
			return $speeds[$i-1];
		}
		else {
			return $speeds[$i];
		}
	}
}

sub divisor {
	# convert dec fraction to 1/x fraction
	# e.g. solve for 1/x = .00095 -> x=1052
	my $t = shift;
	return round_to((1/$t),0);
}

sub exposure {
	# from Covington Appendix A1
	# input f-ratio, film speed, Brightness, filter factor
	# f-ratio expressed as the divisor; if f/16, pass 16
	# B is calculated in sub brightness, without filter factor
	# if no filter factor, pass 1
	# formula t (seconds) = (f**2)/SB
	
	my ($f, $S, $B, $ff) = @_;
	$B = $B/$ff;	#adjusted for filter factor
	return ( ($f**2)/($S*$B) );
}

sub m_from_m2p {
	my ($m2p,$d) = @_;
	return ($m2p - 2.5*( &logbase($pi/4 * $d**2, 10) ));
}

sub B_to_m2p {
	# from Covington Appendix A2
	# m" given brightness
	# brightness and m" are inter-convertible
	my $B = shift;
	return (9.0 - 1.086*log($B));
}

sub brightness {
	# from Covington Appendix A2
	# brightness given calculated m"
	# brightness and m" are inter-convertible
	my $m2p = shift;
	return 2.512**(9.0 - $m2p);
}

sub m2prime {
	# from Covington Appendix A2
	# inputs: total magnitude, apparent diameter in arc-seconds
	# returns magnitude per square arc-second (abbreviated m")
	my ($m, $d) = @_;
	# formula m" = m + 2.5log base 10 (pi/4 * d**2)
	return ($m + 2.5*( &logbase($pi/4 * $d**2, 10) ));
}

sub apparent_diameter {
	# expects distance in AU, radius in km
	# same as Starry Night Pro ephemera
	my ($radius, $distance) = @_;
	$distance = $distance * $au;
	my $diameter = $radius*2;
	my $rad = atan2($diameter,$distance);
	return ( &round_to(($rad*$deg_per_radian*3600),1) );
}

sub logbase { 
	# Example: logbase(243, 3) is 5 -- p. 484 Algorithms
	# log (n, 10) equals log(n)/log(10) - basic algebra
    my ($number, $base) = @_;
    return if $number <= 0 or $base <= 0 or $base == 1;
    return log($number) / log($base);
}

sub round { $_[0] > 0 ? int $_[0] + 0.5 : int $_[0] - 0.5 }

sub round_to {
	my ($raw, $places) = @_;
	my $power = 10**$places;
	return (round($raw *$power)/$power);
}
#============================================================================

1;