Ticket #260: projections.pl

#!/usr/bin/perl

use Math::Trig;

$alpha = pi/4;                  # RA of point of interest
$delta = pi/4;                  # Dec of point of interest
print "alpha,delta: $alpha $delta\n";

$alphaP = 0;                    # RA of tangent point
$deltaP = 0;                    # Dec of tangent point

# Rotate to appropriate system
$theta = asin(sin($delta)*sin($deltaP) + cos($delta)*cos($deltaP)*cos($alpha-$alphaP));
$phi = atan2(-cos($delta)*sin($alpha-$alphaP), sin($delta)*cos($deltaP) - cos($delta)*sin($deltaP)*cos($alpha-$alphaP));
print "phi,theta: $phi $theta\n";

# Projection
$x = cos($theta) * sin($phi) / sin($theta);
$y = - cos($theta) * cos($phi) / sin($theta);
print "x,y: $x $y\n";

# Inverse
$R = sqrt($x**2 + $y**2);
$squiggle = 1/$R;
$phiNew = atan2(-$x,$y);
$thetaNew = atan2($squiggle, 1);
print "phi,theta: $phiNew $thetaNew\n";

# Rotate back to sky
$deltaNew = asin(sin($theta)*sin($deltaP) + cos($theta)*cos($deltaP)*cos($phi));
$alphaNew = $alphaP + atan2(-cos($theta)*sin($phi), sin($theta)*cos($deltaP) - cos($theta)*sin($deltaP)*cos($phi));
print "alpha,delta: $alphaNew $deltaNew\n";