#!/usr/bin/perl
#
#		This script was first written by Alexey Vikhlinin and
#		then modified by Jacob laas (jclaas@gmail.com).
#
# Revision 2016/02/15  jcl
# fixed a bug where the ev photon energy was not interpreted correctly
#
# Revision 2015/09/03  jcl
# added number density to the pressure converter
# fixed a bug the main pressure unit is caught in the regex before the smaller one
#
# Revision 2015/08/20  jcl
# updated "help" argument to allow more options (i.e. -h or --help)
# added photon/energy/mass/distance/temperature/pressure conversions, removed deprecated ones
#
# Revision 2014/07/15  jcl
# refactored some of the code to be more helpful for the lost
#
# Revision 2013/09/13  jcl
# added energy units to more easily convert between wavenum/eV/Eh/kjmol/kcalmol/K
#
# Revision 2013/08/19  jcl
# added option to do "-l <expression>" and force longer print format
#
# Revision 2013/05/13  jcl
# added astronomical coordinate calculator
#
# Revision 2013/05/02  jcl
# converted units to SI and specified units in comments
# fixed physical constants as per "Review of Particle Physics", Berginer et al. (2012) (URL: http://link.aps.org/doi/10.1103/PhysRevD.86.010001)
# added a few examples for offline help
# added wavelength/energy conversions
# added pressures
# added collision frequency & mean free path
# replaced ampere with angstrom
# fixed spelling typos
# added function helpers
#
# /begin{original header by Alexey}
# v 02.06.2007
# Alexey Vikhlinin
# http://hea-www.harvard.edu/~alexey/calc.html
# some code is taken from Math::Trig;
#
# This script is provided 'as is' and without any warranty. For example,
# I'm not paying if you wrap it in a web form and your files get erased:)
# /end{original header by Alexey}
# 

use Term::ReadLine;
use feature 'switch';
no warnings 'experimental::smartmatch';

$log10 = 2.30258509299405;
$pi = 3.14159265358979324; # precision set just a few digits beyond "long" format
#-----------------------------------------------------------------------------
use constant pi    => 3.14159265358979;		# duh
use constant EE    => 2.718281828459045;	# exp(1). Not 'e' because e is charge of electron

use constant Msun  => 1.9885e+30;			# solar mass (kg)
use constant Lsun  => 3.828e+26;			# solar luminosity (W)
use constant Rsun  => 6.9551e+8;			# solar equatorial radius (m)
use constant c     => 299.792458e+6;		# speed of light (m/s)
use constant G     => 6.67384e-11;			# gravitational constant (m^3 / kg / s^2)
use constant e     => 1.602176565e-19;		# charge of electron (C)
use constant h     => 6.62606957e-34; use constant hbar => h/2/pi;   # Planck's constant (J*s)
use constant me    => 9.10938291e-31;		# mass of electron (kg)
use constant mp    => 1.672621777e-27;		# mass of proton (kg)
use constant mu0   => 4 * pi * 1e-7;		# permeability of free space (N / A^2)
use constant ep0   => 1/mu0/c**2;			# permittivity of free space (F / m)
use constant alpha => e**2/2/ep0/h/c;		# fine structure constant
use constant sigmaT=> 6.652458734e-1;		# Thomson cross section=8pi*re^2/3 (barn)
use constant k     => 1.3806488e-23;		# Boltzman constant (J / K)
use constant NA    => 6.02214129e+23;		# Avogadro constant (per mol)
use constant sigma => 5.670373e-8;			# Stefan-Boltzmann constant=pi^2*k^4/60/hbar^3/c^2 (W m^-2 K^-4)
use constant pc    => 3.0856776e+16; use constant kpc => 1000*pc; use constant Mpc => pc*1e6;   # parsec (m)
use constant AU    => 149597870700;			# astronomical unit (m)
use constant ly    => 0.946053e+16;			# light year (m)
#-----------------------------------------------------------------------------
# `second' is used internally
use constant sec   => 1; use constant s => sec;
use constant min   => 60*sec; use constant minutes => min;
use constant hour  => 3600; use constant hr => hour;      # not 'h' because h is Planck's const
use constant day   => 24*hour; use constant yr => 365.242199*day; use constant year=>yr;
use constant MHz   => 1e+6; use constant GHz => 1000*MHz; # frequency
# `meter' is used internally
use constant meter => 1; use constant m => meter; use constant km => 1e3; use constant cm => 1e-2; use constant millimeter => 1e-3; use constant nm => 1e-9;
use constant inch  => 0.0254; use constant ft => 12*inch; use constant feet => ft; use constant yard => 3*ft;
use constant mile  => 5280*ft; use constant mph => mile/hour;
use constant knot  => 0.514444; use constant knots => knot;
use constant A     => 1e-10;                # angstrom (m) <-- was ampere=1
# `kilogram' is used internally
use constant kg    => 1;                    # kilogram
use constant lb    => 0.4535924;            # pound (kg)
use constant oz    => lb/16;
# `Joule' is used internally
use constant J     => 1; use constant mJ => J/1000; # Joule
use constant eV    => 1.602176565e-19; use constant keV => 1000*eV;   # electron volt (J)
use constant wavenum => 8065.54*eV;		# wavenumber, cm^-1 (J)
use constant aJ    => 1e-15;				# attojoule (J)
use constant Eh    => 219474.63*wavenum;	# Hartree (J)
use constant kjmol => 83.5935*wavenum;		# kJ/mol (J)
use constant kcalmol => 349.755*wavenum;	# kcal/mol (J)
use constant K     => 0.695039*wavenum;	# Kelvin (J)
use constant L     => 1e-3;					# liter (in m^3)
use constant amu   => 1.660538921e-27;		# atomic mass unit or dalton (kg)
use constant W     => 1;					# Watt (J/s)
use constant Pa    => 1;					# pascal (kg / m / s^2)
use constant atm   => 1.01325e5; use constant Torr => 133.3224; use constant bar => 1e5; use constant psi => 6.8948e3;
use constant mbar  => bar/1000; use constant mTorr => Torr/1000;
use constant R     => k * NA;				# gas constant (J / K / mol)
use constant Rt    => 62.36367;				# gas constant (m^3 Torr / K / mol)
use constant Jy    => 1e-26* W/meter**2;	# Jansky flux density (per Hz) (10^-23 erg/s/cm^2/Hz)
use constant deg   => pi/180; use constant arcmin => deg/60; use constant arcsec => arcmin/60;
use constant sccm  => 1; use constant scfh => sccm/472; # flow rates
#-----------------------------------------------------------------------------

sub ln {
	log($_[0]);
}
sub lg {
	log($_[0])/$log10;
}
sub fact {
	$s=1;
	for ($i=2;$i<=$_[0];$i++) {
		$s *= $i;
	}
	return $s;
}
sub r2d {
	$_[0]*180.0/$pi
}
sub atan {
	atan2($_[0],1)
}
sub tan {
	my $z = $_[0];
	sin($z)/cos($z)
}
sub acos {
	my $z = $_[0];
	atan2(sqrt(1-$z*$z), $z)
}
sub asin {
	my $z = $_[0];
	atan2($z, sqrt(1-$z*$z))
}
sub C {
	fact($_[0])/(fact($_[1])*fact(($_[0]-$_[1])))
}
sub sind {
	sin($_[0]*pi/180)
}
sub cosd {
	cos($_[0]*pi/180)
}
sub tand {
	tan($_[0]*pi/180)
}
sub asind {
	asin($_[0])*180/pi
}
sub acosd {
	acos($_[0])*180/pi
}
sub atand {
	atan($_[0])*180/pi
}
sub sinh  {
	my $x=$_[0];
	(exp($x)-exp(-$x))/2.0;
}
sub cosh  {
	my $x=$_[0];
	(exp($x)+exp(-$x))/2.0;
}
sub asinh {
	my $x=$_[0];
	log($x+sqrt(1+$x**2))
}
sub acosh {
	my $x=$_[0];
	log($x+sqrt($x**2-1))
}
sub atanh {
	my $x=$_[0];
	0.5*(log(($x+1)/(1-$x)))
}
sub acoth {
	my $x=$_[0];
	atanh(1/$x)
}
sub nint {
	my $x = $_[0]; 
	my $n = int($x);
	if ( $x > 0 ) {
		if ( $x-$n > 0.5) {
			return $n+1;
		}
		else {
			return $n;
		}
	}
	else {
		if ( $n-$x > 0.5) {
			return $n-1;
		}
		else {
			return $n;
		}
	}
}
sub lgamma {  # as per code from numerical recipes
	my $xx = $_[0];
	my $j, $ser, $stp, $tmp, $x, $y;
	my @cof = (0.0, 76.18009172947146, -86.50532032941677,
				24.01409824083091, -1.231739572450155,
				0.1208650973866179e-2, -0.5395239384953e-5);
	my $stp = 2.5066282746310005;
	
	$x = $xx; $y = $x;
	$tmp = $x + 5.5;
	$tmp = ($x+0.5)*log($tmp)-$tmp;
	$ser = 1.000000000190015;
	foreach $j ( 1 .. 6 ) {
		$y+=1.0;
		$ser+=$cof[$j]/$y;
	}
	return $tmp + log($stp*$ser/$x);
}
sub gamma { 
	return exp(&lgamma($_[0]));
}

# begin Jake's custom functions
sub deg2arcsec {
	my $as = $_[0] * 180 / $pi * 3600;
	$as = sprintf("%0.2f",$as);			# reformats to 2 trailing decimal points
	return $as;
}
sub angres {
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[wavelength/frequency],[telescope diameter]\n";
		print "output is: ",
				"[diffraction limited beam size in arcseconds]\n";
		print "     NOTE: wavelength/frequency and diameter must have units, such as:\n",
			"     1.2 mm or 240 GHz, 10.4 m (which should return ~30''; GHz is allowed as a frequency!)\n",
			"     650 nm, 9.25 inch (i.e. 0.5 to 1 as for a large SCT, which you'll never get from your backyard..)\n";
		return NULL;
	}
	$wav = $_[0];
	if ($wav =~ m/(\d+)\w*GHz/) {
		$wav = 300/$1;
	}
	$diam = $_[1];
	$deg = 1.22 * $wav / $diam;
	$arcsec = deg2arcsec($deg);
	return $arcsec;
}
sub photonflux {
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[wavelength in m],[power in W]\n";
		print "output is: ",
				"[flux in photons/s]\n";
		return NULL;
	}
	my $wavelength = $_[0];
	my $energy = h * c / $wavelength;
	my $power = $_[1];
	return $power / $energy;
}
sub pulseflux {
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[pulse power in J],[wavelength in m],[pulse length in s],[beamsize in cm^2]\n";
		print "output is: ",
				"[flux in photons/s]\n";
		return NULL;
	}
	$power = $_[0];
	$wavelength = $_[1];
	$photonenergy = h * c / ($wavelength);
	$time = $_[2];
	$beamsize = $_[3];
	
	return $power / $photonenergy / $time / $beamsize;
}
sub photoyield {
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[abs cross-sec in cm^2],[photon flux density in photons/s/cm^2],[irradiation time in s]\n";
		print "output is: ",
				"[percent yield of photolysis]\n";
		return NULL;
	}
	$sigma = $_[0];
	$fluxdensity = $_[1];
	$time = $_[2];
	
	$photoyield = 1 - exp(-$sigma * $fluxdensity * $time);
	return $photoyield * 100;
}
sub StoT {
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[S in Jy],[freq in GHz],[beamsize in arcsec]\n";
#				"[S in Jy/beam],[freq in GHz],[beamsize in arcsec]\n";
		print "output is: ",
				"[T in K]\n";
		return NULL;
	}
	return $_[0] * 1.22e6 * (1/$_[1])**2 * ($_[2])**(-2);
#	return $_[0] * 13.6 * (300/$_[1])**2 * ($_[2])**(-2);
}
sub TtoS {
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[T in K],[freq in GHz],[beamsize in arcsec]\n";
		print "output is: ",
#				"[S in Jy/beam]\n";
				"[S in Jy]\n";
		return NULL;
	}
#	return $_[0] * ($_[2])**(2) / 13.6 / 300 * $_[1];
	return $_[0] / 1.22e6 * $_[1] * ($_[2])**(2);
}
sub mfp {
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[number density in cm^-3],[particle diam 1 in m]\n";
		print "output is: ",
				"[distance in m]\n";
		print "note the following particle diameters\n",
				"(taken from Hanlon's Vacuum Tech User's Guide, 2nd ed., Appendix B.2)\n",
				"He      2.18 Å\n",
				"H       2.4 Å\n",
				"H2O     2.7 Å\n",
				"H2      2.74 Å\n",
				"CO      3.12 Å\n",
				"C6H6    3.5 Å (ref: http://www.umich.edu/~elements/03chap/html/collision/)\n",
				"O2      3.61 Å\n",
				"År      3.64 Å\n",
				"Åir     3.72 Å   (and average mass of 28.966 Da)\n",
				"N2      3.75 Å\n",
				"CH3OH   4.1 Å\n",
				"CH4     4.14 Å\n",
				"NH3     4.43 Å\n",
				"CO2     4.59 Å\n",
				"(H2O)   4.6 Å\n",
				"C2H4    4.95 Å\n",
				"C2H5OH  5.1 Å\n",
				"C2H6    5.3 Å\n",
				"dust    2e-7 m\n";
		return NULL;
	}
	
	use constant stp => 2.504e19; # sets number density at STP conditions
	
	my $density = $_[0] * 1e6;
	my $diam = $_[1];
	
	return 1/(sqrt(2) * pi * $density * $diam**2);
}
sub vel {
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[mass in kg],[temperature in K]\n"; # add mass and temp..
		print "output is: ",
				"[average kinetic velocity in m/s]\n";
		return NULL;
	}
	my $mass = $_[0];
	my $temp = $_[1];
	return sqrt(8 * k * $temp / pi / $mass);
}
sub collfreq {
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[density1 in cm^-3],[diam1 in m],[mass in kg],[temperature in K]\n";
		print "output is: ",
				"[collision frequency with self in s^-1]\n";
		return NULL;
	}
	use constant stp => 2.504e19; # sets number density at STP conditions
	
	my $density = $_[0];
	my $diam = $_[1];
	my $mass = $_[2];			# reduced mass
	my $temp = $_[3];
	
	my $mfp = &mfp($density,$diam);
	my $velocity = &vel($mass,$temp);
	return $velocity / $mfp;
}
sub spheremass {
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[radius in m],[density in g/cm^3]\n";
		print "output is: ",
				"[mass in kg]\n";
		return NULL;
	}
	my $rad = $_[0];
	my $vol = 4/3 * pi * $rad**3;
	my $rho = $_[1];
	my $mass = $vol * $rho * 1e3;
	return $mass;
}
sub collfreqism {
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[total density in cm^-3],[fractional abund],[diam in m],[mass in kg],[temperature in K],[dust or H2]\n";
		print "optional arguments are: ",
				"[dust temp],[fraction of H2 wrt total H]\n";
		print "output is: ",
				"[collision frequency with dust in s^-1]\n";
		return NULL;
	}
	use constant dustdensity => 2.5;		# mass density of dust in g/cm^3 (ranges from 2-3)
	use constant dustradius => 1e-7;		# radius of dust grain
	my $massdust = &spheremass(dustradius,dustdensity);	# mass of dust in kg
	
	my ($fracdust);
	if ($_[7]) {
		$fracdust = ($_[7] * 4 + (1 - 2*$_[7])) * amu / $massdust * 0.01; # should be approx 10^-12
	} else {
		$fracdust = 2 * amu / $massdust * 0.01;	# assumes all hydrogen is as H2
	}
	
	my $densitytot = $_[0];
	my $densitydust = $densitytot * $fracdust;
	my $fracgas = $_[1];
	my $densitygas = $densitytot * $fracgas;
	my $diamgas = $_[2];
	my $massgas = $_[3];
	my $tempgas = $_[4];
	my ($tempdust);
	if ($_[6]) {
		$tempdust = $_[6];
	} else {
		$tempdust = $tempgas;
	}
	if ($_[7]) {
		$densityH = $densitytot * (1-2*$_[7]);
		$densityH2 = $densitytot * $_[7];
	} else {
		$densityH = 0;
		$densityH2 = $densitytot * 0.5;
	}
	
	my $veldust = &vel($massdust,$tempdust);
	my $velgas = &vel($massgas,$tempgas);
	my $velH = &vel(amu,$tempgas);
	my $velH2 = &vel(3.32108e-27,$tempgas);
	my $velrel = sqrt($velH2**2 + $velgas**2);
	print	"the following velocities were determined:\n",
			"          dust - $veldust m/s\n",
			"           gas - $velgas m/s\n",
			"relative to H2 - $velrel m/s\n";
	
	my $zmolecwithself = pi * $diamgas**2 * sqrt(2)*$velgas * $densitygas;
	#my $bmax = ($diamgas + 2.74 * A) / 2;
	my $zmolecwithH = pi * (($diamgas + 2.5 * A)/2)**2 * $velrel * $densityH;
	my $zmolecwithH2 = pi * (($diamgas + 2.74 * A)/2)**2 * $velrel * $densityH2;
	my $zmolecwithdust = pi * dustradius**2 * $velgas * $densitydust;
	my $ztotpervol = $zmolecwithdust * $densitygas / 2;
	
	printf "coll with self is:   %.5e\n",$zmolecwithself;
	printf "coll with H is:      %.5e\n",$zmolecwithH;
	printf "coll with H2 is:     %.5e\n",$zmolecwithH2;
	printf "coll with dust is:   %.5e\n",$zmolecwithdust;
	
	if ($_[5] =~ /dust/) {
		return $zmolecwithdust;
	} else {
		return $zmolecwithH2;
	}
}
sub BtoJ { # uses process outlined at http://www.stargazing.net/kepler/b1950.html
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[RA in hr,min,sec],[DEC in deg,min,sec]\n";
		print "output is: ",
				"[RA in hr:min:sec],[DEC in deg:min:sec]\n";
		return NULL;
	}
	my $tempdeg = $_[0] * 15;
	my $rad = ($_[0] + $_[1]/60 + $_[2]/3600) * 15;
	my $decsign = 1;
	if ($_[3] < 0) {
		$decsign = -1;
	}
	my $decd = $_[3] + $_[4]/60*$decsign + $_[5]/3600*$decsign;
	
	# warn about limits of accuracy near poles
	if (abs($decd) < 1 or abs($decd) > 88) {
		print "WARNING: This is way too close to a pole -- coordinates may be off by up to 2 degrees\!\n\n";
		print "Use should use SIMBAD (http://simbad.u-strasbg.fr/simbad/sim-fcoo)..\n";
	} elsif (abs($decd) < 2 or abs($decd) > 85) {
		print "WARNING: This is pretty close to a pole -- coordinates may be off by up to 10\'\'\!\n\n";
	} elsif (abs($decd) < 3 or abs($decd) > 80) {
		print "DEC is fairly close to a pole -- coordinates may be off by up to 5\'\'\!\n\n";
	} else {}
	print "SIMBAD query is: http://simbad.u-strasbg.fr/simbad/sim-coo?CooDefinedFrames=none&CooEpoch=1950&Coord="."$rad"."%20"."$decd"."&submit=submit%20query&Radius.unit=arcmin&CooEqui=1950&CooFrame=FK5&Radius=2 \n\n";
	
	# convert to rectangular coordinates
	my $x = cosd($rad) * cosd($decd);
	my $y = sind($rad) * cosd($decd);
	my $z = sind($decd);
	my @orig = ($x, $y, $z);
	my @new = (0,0,0);
	
	# for EQ 1950, conversion is simple using the transfer matrix here (see S3.591 of Explanatory Supplement 2005 for ):
	my @matrix = (
		[0.9999256782, -0.0111820611, -0.0048579477],
		[0.0111820610,  0.9999374784, -0.0000271765],
		[0.0048579479, -0.0000271474,  0.9999881997]
	);
	for (my $i=0;$i<3;$i++) {
		for (my $j=0;$j<3;$j++) {
			$new[$i] += $orig[$j] * $matrix[$i][$j];
		}
	}
	
	# convert back to new ra/dec
	#$newrad = $rad + 0.640265 + 0.278369 * sind($rad) * tan($decd); # lower precision version
	my $newrad = atan($new[1] / $new[0]) * 180 / pi;
	if ($new[0] < 0) {
		$newrad += 180;
	} elsif ($new[1] < 0 && $new[0] > 0) {
		$newrad += 360;
	}
	#my $newdecd = $decd + 0.278369 * cosd($rad); # lower precision
	my $newdecd = asin($new[2]) * 180 / pi; # rect coordinate conversion
	
	my $newrahr = int($newrad/15);
	my $min = ($newrad/15 - $newrahr) * 60;
	my $newramin = int($min);
	my $newrasec = ($newrad/15 - $newrahr - $newramin/60) * 3600;
	$newrasec = sprintf("%0.2f",$newrasec); # cut back on precision
	my $newdec = int($newdecd);
	my $min = ($newdecd - $newdec) * 60 * $decsign;
	my $newdecmin = int($min);
	my $newdecsec = ($newdecd - $newdec - $newdecmin/60*$decsign) * 3600 * $decsign;
	$newdecsec = sprintf("%0.2f",$newdecsec); # cut back on precision
	
	print "J2000 RA/DEC -- $newrahr:$newramin:$newrasec $newdec:$newdecmin:$newdecsec\n";
	return 1;
}
sub Jday { # note J2000.0 is defined as 2451545.0
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[year],[numerical month],[numerical day]\n";
		print "output is: ",
				"[Julian day]\n";
		return NULL;
	}
	my $a = int($_[0]/100);
	my $b = int($a/4);
	my $c = int(2 - $a + $b);
	my $e = int(365.25 * ($_[0]+4716));
	my $f = int(30.6001 * ($_[1]+1));
	my $Jday = $c + $_[2] + $e + $f - 1524.5;
	return $Jday;
}
sub degtodecimal { # works for deg:min:sec
	my $sign = 1;
	if ($_[0] < 0) {
		$sign = -1;
	}
	return $_[0] + ($_[1]/60 + $_[2]/3600)*$sign;
}
sub Jtohorizon {
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[RA in hr,min,sec],4,[DEC in deg,min,sec],7,[LAT in deg,min,sec],11,[LONG in deg,min,sec],15,[date in year,mon,day],19,[time in hr,min,sec]\n",
				"RA/DEC should be in J2000.0, date/time is local, and time should be in long (24hr) format; note the placeholders between groups\n";
		print "output is: ",
				"[altitude in deg],[azimuth in deg]\n";
		return NULL;
	}
	
	# convert to deg
	my $tempdeg = $_[0] * 15; # converts RA's hr to deg first...
	my $rad = &degtodecimal($tempdeg,$_[1],$_[2]);
	my $decd = &degtodecimal($_[4],$_[5],$_[6]);
	my $latd = &degtodecimal($_[8],$_[9],$_[10]);
	my $longd = &degtodecimal($_[12],$_[13],$_[14]);
	
	my $lh = $_[20] + $_[21]/60 + $_[22]/3600; # local hour
	my $UT = $lh + $longd/15;
	if ($UT > 24) {
		$UT -= 24;
	}
	my $JD = &Jday($_[16], $_[17], $_[18]) + $lh/24;
	my $sincedays = $JD - 2451545; # Julian days since J2000.0
	print "$sincedays days since J2000.0\n";
	my $lst = 100.46 + 0.985647 * $sincedays;# + $longd + 15 * $UT;
	if ($lst < 360) {
		while ($lst < 0) {$lst +=360;}
	} else {
		while ($lst > 360) {$lst -= 360;}
	}
	print "LST deg is: $lst\n";
	
	my $ha = $lst - $rad;
	if ($ha < 0) {$ha += 360;}
	print "hour angle is $ha\n";
	
	$ha = 54.382617;
	$decd = 36.466667;
	$latd = 52.5;
	my $sinalt = sind($decd) * sind($latd) + cosd($decd) * cosd($latd) * cosd($ha);
	my $alt = asin($sinalt);
	my $cosaz = sind($decd) - sin($alt) * sind($latd) / cos($alt) / cosd($latd);
	my $az = acos($cosaz);
	
	$alt = $alt * 180/pi;
	$az = $az * 180/pi;
	print "alt/alt are $alt and $az\n";
	
	return 1;
}
sub vishorizon {
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[height of observer in ft]\n";
		print "output is: ",
				"[distance to horizon in miles]\n";
		return NULL;
	}
	$height = $_[0]/5286; # to miles
	$earthrad = 3958.76;# in miles
	$earthdiam = 2*$earthrad;# in km
	return sqrt($height*($earthdiam+$height));
}
sub visobject {
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[height of observer in ft],[height of object in ft]\n";
		print "output is: ",
				"[maximum viewable distance in miles]\n";
		return NULL;
	}
	$d_v = vishorizon($_[0]);
	$d_obj = vishorizon($_[1]);
	$maxdistance = $d_v + $d_obj;
	print "Add ~8% if you want to consider visible refraction under standard conditions!\n";
	print "Add a total of 15% if you want to consider radar (30-300 mm)...\n";
	return $maxdistance;
}
sub halflife {
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[half-life],[time which has passed]\n";
		print "output is: ",
				"[remaining fraction (N/N0) of material]\n";
		return NULL;
	}
	return 2**(-$_[1]/$_[0]);
}
sub drake {
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[R],[fp],[ne],[fl],[fi],[fc],[L], where:\n",
				"R is the average rate of star formation in the galaxy\n",
				"fp is the fraction of suitable suns with planetary systems\n",
				"ne is the mean number of planets that are located within the zone where water can exist as a liquid\n",
				"fl is the fraction of such planets on which life actually originates\n",
				"fi represents the fraction of such planets on which some form of intelligence arises\n",
				"fc is the fraction of such intelligent species that develop the ability and desire to communicate with other civilizations, and\n",
				"L is the mean lifetime (in years) of a communicative civilization\n";
		print "Originally: (1, 0.2-0.5, 1-5, 1, 1, 0.1-0.2, 1e3-1e8)\n",
				"current estimates: (7, 1.6, 0.3, 0.13, 1, 1, 5e5)\n";
		print "output is: ",
				"[the no. of active, communicative civ in the MW]\n";
		return NULL;
	}
	return $_[0]*$_[1]*$_[2]*$_[3]*$_[4]*$_[5]*$_[6];
}
sub temperature {
	# see https://en.wikipedia.org/wiki/Conversion_of_units_of_temperature#Kelvin
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[value],[unit],[desired unit]\n";
		print "output is: ",
				"either a) a list of conversions + return NULL, or b) only the new unit\n";
		print "the following units are available (not case-sensitive):\n",
				"K - Kelvin (SI)\n",
				"C - Celsius\n",
				"F - Fahrenheit\n",
				"Ra - Rankine\n",
				"D - Delisle\n",
				"N - Newton\n",
				"Re - Réaumur\n",
				"Ro - Rømer\n";
		return NULL;
	}
	# define preferred unit (e.g. SI)
	my $kelvin;
	# test for cases and convert to primary unit
	given (lc $_[1]) {
		when(/k/) { $kelvin = $_[0] }
		when(/c/) { $kelvin = $_[0] + 273.15 }
		when(/f/) { $kelvin = ($_[0] + 459.67) * 5/9 }
		when(/ra/) { $kelvin = $_[0] * 5/9 }
		when(/d/) { $kelvin = 373.15 - $_[0] * 2/3 }
		when(/n/) { $kelvin = $_[0] * 100/33 + 273.15 }
		when(/re/) { $kelvin = $_[0] * 5/4 + 273.15 }
		when(/ro/) { $kelvin = ($_[0] - 7.5) * 40/21 + 273.15 }
		default {
			print "you did not specify a unit!\n";
			temperature("help");
			return NULL;
		}
	}
	# convert from preferred unit
	my $celsius = $kelvin - 273.15;
	my $fahrenheit = $kelvin * 9/5 - 459.67;
	my $rankine = $kelvin * 9/5;
	my $delisle = (373.15 - $kelvin) * 3/2;
	my $newton = ($kelvin - 273.15) * 33/100;
	my $reaumer = ($kelvin - 273.15) * 4/5;
	my $romer = ($kelvin - 273.15) * 21/40 + 7.5;
	# if a specific unit is desired, return it
	if (lc $_[2]) {
		given (lc $_[2]) {
			when(/k/) { return $kelvin }
			when(/c/) { return $celsius }
			when(/f/) { return $fahrenheit }
			when(/ra/) { return $rankine }
			when(/d/) { return $delisle }
			when(/n/) { return $newton }
			when(/re/) { return $reaumer }
			when(/ro/) { return $romer }
			default {
				print "the requested unit could not be correctly interpreted\n";
				temperature("help");
				return NULL;
			}
		}
	}
	# ..otherwise print all the units out
	print "$_[0] $_[1] is:\n",
			"\t$kelvin Kelvin\n",
			"\t$celsius °C\n",
			"\t$fahrenheit °F\n",
			"\t$rankine °R\n",
			"\t$delisle °De\n",
			"\t$newton °N\n",
			"\t$reaumer °Ré\n",
			"\t$romer °Rø\n";
	return NULL;
}
sub distance {
	# see http://www.convert-me.com/en/convert/length/meter.html
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[value],[unit],[desired unit]\n";
		print "output is: ",
				"either a) a list of conversions + return NULL, or b) only the new unit\n";
		print "the following units are available (not case-sensitive):\n",
				"m - meter (SI)\n",
				"i - inch\n",
				"ft - foot\n",
				"yd - yard\n",
				"mile - mile\n",
				"AU - astronomical unit\n",
				"ly - light-year\n",
				"pc - parsec\n";
		return NULL;
	}
	# define preferred unit (e.g. SI)
	my $meter;
	# test for cases and convert to primary unit
	given (lc $_[1]) {
		when(/m/) { $meter = $_[0] }
		when(/i/) { $meter = $_[0] / 39.370079 }
		when(/ft/) { $meter = $_[0] / 3.2808399 }
		when(/yd/) { $meter = $_[0] / 1.0936133 }
		when(/mile/) { $meter = $_[0] / 0.00062137119 }
		when(/au/) { $meter = $_[0] * 149597870700 }
		when(/ly/) { $meter = $_[0] * 9460730472580800 }
		when(/pc/) { $meter = $_[0] * 3.085677581e16 }
		default {
			print "you did not specify a unit!\n";
			distance("help");
			return NULL;
		}
	}
	# convert from preferred unit
	my $inch = $meter * 39.370079;
	my $foot = $meter * 3.2808399;
	my $yard = $meter * 1.0936133;
	my $mile = $meter * 0.00062137119;
	my $au = $meter / 149597870700;
	my $ly = $meter / 9460730472580800;
	my $pc = $meter / 3.085677581e16;
	# if a specific unit is desired, return it
	if ($_[2]) {
		given (lc $_[2]) {
			when(/m/) { return $meter }
			when(/i/) { return $inch }
			when(/ft/) { return $foot }
			when(/yd/) { return $yard }
			when(/mile/) { return $mile }
			when(/au/) { return $au }
			when(/ly/) { return $ly }
			when(/pc/) { return $pc }
			default {
				print "the requested unit could not be correctly interpreted\n";
				distance("help");
				return NULL;
			}
		}
	}
	# ..otherwise print all the units out
	print "$_[0] $_[1] is:\n",
			"\t$meter m\n",
			"\t$inch in.\n",
			"\t$foot ft\n",
			"\t$yard yd\n",
			"\t$mile mi\n",
			"\t$au AU\n",
			"\t$ly ly\n",
			"\t$pc pc\n";
	return NULL;
}
sub mass {
	# see http://www.convert-me.com/en/convert/weight/kilogram.html
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[value],[unit],[desired unit]\n";
		print "output is: ",
				"either a) a list of conversions + return NULL, or b) only the new unit\n";
		print "the following units are available (not case-sensitive):\n",
				"kg - kilogram (SI)\n",
				"lb - pound\n",
				"oz - ounce\n";
		return NULL;
	}
	# define preferred unit (e.g. SI)
	my $kg;
	# test for cases and convert to primary unit
	given (lc $_[1]) {
		when(/kg/) { $kg = $_[0] }
		when(/lb/) { $kg = $_[0] / 2.2046226 }
		when(/oz/) { $kg = $_[0] / 35.273962 }
		default {
			print "you did not specify a unit!\n";
			mass("help");
			return NULL;
		}
	}
	# convert from preferred unit
	my $lb = $kg * 2.2046226;
	my $oz = $kg * 35.273962;
	# if a specific unit is desired, return it
	if ($_[2]) {
		given (lc $_[2]) {
			when(/kg/) { return $kg }
			when(/lb/) { return $lb }
			when(/oz/) { return $oz }
			default {
				print "the requested unit could not be correctly interpreted\n";
				mass("help");
				return NULL;
			}
		}
	}
	# ..otherwise print all the units out
	print "$_[0] $_[1] is:\n",
			"\t$kg kg\n",
			"\t$lb lbs\n",
			"\t$oz oz.\n";
	return NULL;
}
sub pressure {
	# see http://www.convert-me.com/en/convert/pressure/pascal.html
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[value],[unit],[desired unit]\n";
		print "output is: ",
				"either a) a list of conversions + return NULL, or b) only the new unit\n";
		print "the following units are available (not case-sensitive):\n",
				"Pa - Pascal (SI)\n",
				"bar - bar\n",
				"atm - physical atmosphere\n",
				"Torr - Torr (aka mmHg)\n",
				"mTorr - milliTorr\n",
				"psi - lbs per sq. inch\n",
				"pcc - per cm^3 (at 20°C)\n";
		return NULL;
	}
	# define preferred unit (e.g. SI)
	my $pascal;
	# test for cases and convert to primary unit
	given (lc $_[1]) {
		when(/pa/) { $pascal = $_[0] }
		when(/mbar/) { $pascal = $_[0] * 100 }
		when(/bar/) { $pascal = $_[0] * 1e5 }
		when(/atm/) { $pascal = $_[0] * 101325 }
		when(/mtorr/) { $pascal = $_[0] * 0.133322 }
		when(/torr/) { $pascal = $_[0] * 133.322 }
		when(/psi/) { $pascal = $_[0] / 0.00014503774 }
		when(/pcc/) { $pascal = $_[0] * k * 293.15 * 1e-6 }
		default {
			print "you did not specify a unit!\n";
			pressure("help");
			return NULL;
		}
	}
	# convert from preferred unit
	my $bar = $pascal / 1e5;
	my $mbar = $pascal / 100;
	my $atm = $pascal / 10120;
	my $torr = $pascal / 133.32;
	my $mtorr = $pascal / 0.13332;
	my $psi = $pascal / 6894.6;
	my $pcc20C = $pascal / k / 293.15 / 1e6;
	my $pcc10K = $pascal / k / 10 / 1e6;
	# if a specific unit is desired, return it
	if ($_[2]) {
		given (lc $_[2]) {
			when(/pa/) { return $pascal }
			when(/mbar/) { return $mbar }
			when(/bar/) { return $bar }
			when(/atm/) { return $atm }
			when(/mtorr/) { return $mtorr }
			when(/torr/) { return $torr }
			when(/psi/) { return $psi }
			when(/pcc/) { return $pcc20C }
			default {
				print "the requested unit could not be correctly interpreted\n";
				pressure("help");
				return NULL;
			}
		}
	}
	# ..otherwise print all the units out
	print "$_[0] $_[1] is:\n",
			"\t$pascal Pa\n",
			"\t$bar bar\n",
			"\t$mbar mbar\n",
			"\t$atm atm\n",
			"\t$torr Torr\n",
			"\t$mtorr mTorr\n",
			"\t$psi psi\n",
			"\t$pcc20C cm^-3 (at 20°C)\n",
			"\t$pcc10K cm^-3 (at 10K)\n";
	return NULL;
}
sub energy {
	# see http://www.hbcpnetbase.com/
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[value],[unit],[desired unit]\n";
		print "output is: ",
				"either a) a list of conversions + return NULL, or b) only the new unit\n";
		print "the following units are available (not case-sensitive):\n",
				"J - Joule (SI)\n",
				"wavenum - wavenumber (inverse cm)\n",
				"K - Kelvin\n",
				"eV - electronvolts\n",
				"Eh - Hartree\n",
				"kjmol - kJ/mol\n",
				"kcalmol - kcal/mol\n";
		return NULL;
	}
	# define preferred unit (e.g. SI)
	my $joule;
	# test for cases and convert to primary unit
	given (lc $_[1]) {
		when(/j/) { $joule = $_[0] }
		when(/wavenum/) { $joule = $_[0] * 1.986447e-5 * 1e18 }
		when(/ev/) { $joule = $_[0] * 0.1602177 * 1e18 }
		when(/eh/) { $joule = $_[0] * 4.359748 * 1e18 }
		when(/kjmol/) { $joule = $_[0] * 1.660540e-3 * 1e18 }
		when(/kcalmol/) { $joule = $_[0] * 6.9477e-3 * 1e18 }
		when(/k/) { $joule = $_[0] * 1.380658e-5 * 1e18 }
		default {
			print "you did not specify a unit!";
			energy("help");
			return NULL;
		}
	}
	# convert from preferred unit
	my $wavenum = $joule * 50341.1 * 1e-18;
	my $ev = $joule * 6.241506 * 1e-18;
	my $eh = $joule * 0.2293710 * 1e-18;
	my $kjmol = $joule * 602.2137 * 1e-18;
	my $kcalmol = $joule * 143.9325 * 1e-18;
	my $kelvin = $joule * 7.24292e4 * 1e-18;
	# if a specific unit is desired, return it
	if ($_[2]) {
		given (lc $_[2]) {
			when(/j/) { return $joule }
			when(/wavenum/) { return $wavenum }
			when(/ev/) { return $ev }
			when(/eh/) { return $eh }
			when(/kjmol/) { return $kjmol }
			when(/kcalmol/) { return $kcalmol }
			when(/k/) { return $kelvin }
			default {
				print "the requested unit could not be correctly interpreted\n";
				energy("help");
				return NULL;
			}
		}
	}
	# ..otherwise print all the units out
	print "$_[0] $_[1] is:\n",
			"\t$joule Joules\n",
			"\t$wavenum cm-1\n",
			"\t$ev eV\n",
			"\t$eh Hartrees\n",
			"\t$kjmol kj/mol\n",
			"\t$kcalmol kcal/mol\n",
			"\t$kelvin K\n";
	return NULL;
}
sub photon {
	# see https://en.wikipedia.org/wiki/Conversion_of_units_of_temperature#Kelvin
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[value],[unit],[desired unit]\n";
		print "output is: ",
				"either a) a list of conversions + return NULL, or b) only the new unit\n";
		print "the following units are available (not case-sensitive):\n",
				"Hz - Hertz (SI)\n",
				"these frequencies: MHz, GHz, THz\n",
				"these wavelengths: nm, um, mm\n",
				"wav - wavenumber (inverse cm)\n";
				"ev - electronvolts\n";
		return NULL;
	}
	# define preferred unit (e.g. SI)
	my $hertz;
	# test for cases and convert to primary unit
	given (lc $_[1]) {
		when(/mhz/) { $hertz = $_[0] * 1e6 }
		when(/ghz/) { $hertz = $_[0] * 1e9 }
		when(/thz/) { $hertz = $_[0] * 1e12 }
		when(/nm/) { $hertz = c / $_[0] * 1e9 }
		when(/um/) { $hertz = c / $_[0] * 1e6 }
		when(/mm/) { $hertz = c / $_[0] * 1e3 }
		when(/wav/) { $hertz = $_[0] * c * 100 }
		when(/ev/) { $hertz = $_[0] * 2.417988e14 }
		when(/hz/) { $hertz = $_[0] }	# because the regex is greedy
		default {
			print "you did not specify a unit!\n";
			photon("help");
			return NULL;
		}
	}
	# convert from preferred unit
	my $mhz = $hertz / 1e6;
	my $ghz = $hertz / 1e9 ;
	my $thz = $hertz / 1e12 ;
	my $nm = c / $hertz * 1e9;
	my $um = c / $hertz * 1e6;
	my $mm = c / $hertz * 1e3;
	my $wav = $hertz / c / 100;
	my $ev = $hertz * 4.13558e-15;
	# if a specific unit is desired, return it
	if ($_[2]) {
		given (lc $_[2]) {
			when(/mhz/) { return $mhz }
			when(/ghz/) { return $ghz }
			when(/thz/) { return $thz }
			when(/nm/) { return $nm }
			when(/um/) { return $um }
			when(/mm/) { return $mm }
			when(/wav/) { return $wav }
			when(/ev/) { return $ev }
			when(/hz/) { return $hertz }
			default {
				print "the requested unit could not be correctly interpreted\n";
				photon("help");
				return NULL;
			}
		}
	}
	# ..otherwise print all the units out
	print "$_[0] $_[1] is:\n",
			"\t$hertz Hz\n",
			"\t$mhz MHz\n",
			"\t$ghz GHz\n",
			"\t$thz THz\n",
			"\t$nm nm\n",
			"\t$um μm\n",
			"\t$mm mm\n",
			"\t$wav cm-1\n",
			"\t$ev eV\n";
	return NULL;
}
sub placeholder {
	if ($_[0] =~ /help/) {
		print "arguments should be: ",
				"[insert something helpful here]\n";
		print "optional arguments are: ",
				                      "[blah1]",
			  "                        [blah2]\n";
		print "output is: ",
				"[insert something helpful here]\n";
		return NULL;
	}
	return NULL;
}
# end Jake's custom functions

%SHELLCOM=();
do $ENV{"CALCRC"} if ($ENV{"CALCRC"}); # personal settings

sub doeval {
	# remove "calc " from the beginning of the line
	s/^\s*calc\s//;
	$origexpr = $_;
	
	# Jan 19 2007: the shell is likely to put spaces around ( ); 
	# for us this is unnecessary and indeed, harmful 
	# (breaks the 45mp -> (45*mp) subst.) So remove those spaces
	s/(^| )([()])( |$)/$2/g;
	
	# Jan 19 2007: I've seen in somebody's examples calc 6^2 meaning 6**2;
	# this is wrong! Perl does something else to 6^2, so
	s/\^/**/g;
	
	# Used to be: s/(^|[^\w\+-])....
	s/(^|[^\w])([\+-]?(?:\d+\.?\d*|\.\d+))([eE])([\+-]?\d+(?=$|[^\d\.]))/$1$2\371$4/g;
	# 1.5e5 -> 1.5�5 so that scientific notations don't mess up with the charge of electron
	s/(^\'|\'$)//g; # Trim single quotes in the beg/end of string
	
	s/(^|[^a-zA-Z])s($|[^a-zA-Z])/$1sec$2/g;  # s -> sec
	s/(^|[^a-zA-Z])mm($|[^a-zA-Z])/$1millimeter$2/g;# mm -> millimeter
	s/(^|[^a-zA-Z])m($|[^a-zA-Z])/$1meter$2/g;# m -> meter
	s/(\W)([a-zA-Z]\w*)!/$1fact($2)/g;  # $na_3! -> fact($n)
	
	# Mon Feb  5 09:54:58 2007:  also, 10.(me+mp) -> 10.*(me*mp)
	s/(\d|\.)\s*\(/$1*\(/g;            # 10(me+mp) -> 10*(me+mp)
	# #  s/(^|[^\w\371\+-])([\+-]?(?:\d+\.?\d*|\.\d+)(?:\371[\+-]?\d+)?)([a-zA-Z]+)(?=$|\W)/$1\($2*$3\)/g;              # 45mp -> (45*mp) (to be able to 6me/3mp)
	#  s/(^|[^\w\371])((?:\d+\.?\d*|\.\d+)(?:\371[\+-]?\d+)?)([a-zA-Z]+)(?=[\$\/\*-\+\s]|$)/$1\($2*$3\)$4/g;              # 45mp -> (45*mp) (to be able to 6me/3mp)
	# Modified prev line on Jan 19 2007; was:
	# s/(^|[^\w\371])((?:\d+\.?\d*|\.\d+)(?:\371[\+-]?\d+)?)([a-zA-Z]+)(?=$|\W)/$1\($2*$3\)/g;              # 45mp -> (45*mp) (to be able to 6me/3mp)
	# otherwise asind(0.5sqrt(2)) was broken
	
	s/([\d\.\)])\s*([a-zA-Z])/$1*$2/g; # 45a -> 45*a; )a-> )*a
	s/(\w)\s+([a-zA-Z\d])/$1*$2/g;  # Msun c**2 -> Msun*c**2
	s/(^|\W)([abdfgijln-rt-z])(?=$|\W)/$1\$$2/g;	# single letter lower-case variables
	s{/([abdfgijln-r-t-z])(?=$|\W)}{/\$$1}g;		# preceeded and followed by a 
													# non-word character are replaced
													# with $variable; second substitution
													# is needed in cases /n, because $1
													# works incorrectly
	s/(\d+)!/fact($1)/g; # 10! -> fact(10)
	# factorial after the numbers is done last because otherwise the typo
	# "calc 3.5!" results (silently) in 3.*fact(5)
	s/(^|[^\w\.])0+(\d)/$1$2/g; # trim leading zeros in numbers
	
	s/\:/\,/g; # convert colons to commas (useful for astro coordinates)
	
	s/\371/e/g; # 1.5�5 -> 1.5e5 THIS MUST BE THE LAST TRANSFORMATION
	
	if ($debug) {print STDERR "Debug: ",$_," format: ",$format,"\n";}
	if ( !defined($_result_ = eval "$_") ) {
		print STDERR "wrong expr: $origexpr -> $_\n";
	}
	
	if ( $commandline || (! /=/) ) { # ignore assignments in the stdin mode
		if ($format =~ /^%l/ ) {print $_result_,"\n";}
		else {printf "$format\n", $_result_;}
	}
}

# ENV VARIABLES
sub check_env {
	if ($ENV{"CALCFORMAT"}) {$format = $ENV{"CALCFORMAT"};}
	else {$format = "%g"}
	if ($ENV{"CALCDEBUG"}) {$debug=1}
	else {$debug=0}
}

$_ = "@ARGV";
&check_env;

if ( /^(-h|-{0,2}help)$/ ) {
	&printconst;
	exit 0;
}

if ( /^(-functions)$/ ) {
	exit 0;
}

if ( /\w/ ) {
	$commandline = 1;
	if (/^\s*(-l)(.*)$/) {
		print "printing long format..\n";
		$format = "%l";
		$_ = $2;
		&doeval;
	} else {
		&doeval;
	}
}
else {
	$term = new Term::ReadLine 'Perl calc';
	
	#  print STDERR $term->ReadLine,"\n";
	
	$prompt = "calc: ";
	$OUT = $term->OUT || STDOUT;
	while ( defined ($_ = $term->readline($prompt)) ) {
		if (/^\s*history\s*$/ ) {
			@govgov=$term->GetHistory;
			foreach $gov ( @govgov ) {
				if (! ($gov =~ /^\s*history\s*$/) ) {
					print $gov,"\n";
				}
			}
			next;
		}
		if (/^\s*calc\s*$/ ) {exec $0;}
		if (/^\s*setenv CALC/) {
			@gov = split;
			$ENV{$gov[1]}=$gov[2];
			&check_env;
			next;
		}
		if (/^\s*unsetenv CALC/) {
			@gov = split;
			$ENV{$gov[1]}=0;
			&check_env;
			next;
		}
		if (/^\s*(-h|-{0,2}help)\s*$/) {
			&printconst;
			next;
		}
		if (/^\s*(-l)(.*)$/) {
			print "using long format..\n";
			$format = "%l";
			$_ = $2;
			$commandline = 0;
			&doeval;
			$format = "%g";
			next;
		}
		@gov=split;
		if ($SHELLCOM{$gov[0]}) {
			$gov=shift(@gov);
			system($SHELLCOM{$gov},@gov);
			next;
		}
		$commandline = 0;
		&doeval;
	}
#	while ( <STDIN> ) {
#		chop;
#		$commandline = 0;
#		&doeval;
#	}
}

sub printconst {
	$scriptname = $0;
	print "Available constants (all units are in the SI system):\n";
	open (SCR,$scriptname) or die "Can't open $scriptname";
	while (<SCR>) {
		if (/^use constant |^\#------/) {
			s/use constant //g;
			s/=>/=/g;
			print;
		}
	}
	while (<DATA>) { print; }
}

__END__
Available functions:
-------------------
exp, log, abs, sqrt, sin, cos, tan, asin, acos, atan, atan2(x,y)=atan(x/y)
sind, cosd, tand, asind, acosd, atand --- degree based trig. functions
sinh(x), cosh(x), asinh(x), acosh(x), atanh(x), acoth(x) --- hyperbolic trig.
ln(x), lg(x) --- natural and base-10 logarithms
gamma(x), lgamma(x) --- Gamma and log(Gamma) functions
r2d(x)  = x*180/pi;  or you can use the 'deg' constant
fact(n) = n!; C(n,k) = n!/(k!(n-k)!)
int, nint --- integer portion and nearest integer
photonflux, pulseflux, StoT, TtoS --- deals with photon flux, power and brightness
photoyield --- determines time-dependent photolysis yield (one-time use?)
mfp, vel, collfreq, collfreqism --- gas kinetics
angres --- returns the diffraction limited beam size of a telescope
BtoJ --- converts B1950.0 to J2000.0
Jtohorizon --- convert J2000.0 coordinate to alt/az in absolute degrees
Jday --- determines the Julian day of the year/month/day
vishorizon, visobject --- visible distances to horizon and tall objects
halflife --- calculates the remaining fraction
drake --- the drake equation!
temperature, pressure --- temperature and pressure conversions
mass, distance --- mass and distance conversions
energy, photon --- convert energy units and photon energies/wavelengths

Examples of usage can be found on http://hea-www.harvard.edu/~alexey/calc-examples.html
Some examples include:
% calc 2*(2/8) + 3**2			# simple arithmetics
% calc sin(30deg)				# sin and asin in degrees
% calc 412 pc/ly				# distance to Orion Nebula in light years
% calc temperature(23,'C')		# convert 23 degrees C
% calc photon(185,'nm')			# see energetics of 185 nm photon
% calc pressure(400,'Torr')		# convert 400 Torr to other pressure units

TO ADD:
- http://en.wikipedia.org/wiki/Doomsday_rule
- pH from buffer
- radial velocity (http://ircamera.as.arizona.edu/astr_250/Lectures/Lec_22sml.htm)
- distances of astronomical coordinates
- other astro calcs - http://129.79.46.40/~foxd/cdrom/dos/manual10.htm