#!/usr/local/bin/perl

use Math::Trig;

$DEG_IN_RAD = 57.2957795130823;

# l2,b2
$alphaP = (90 + 192.85948) / $DEG_IN_RAD;
$deltaP = (90 - 27.12825) / $DEG_IN_RAD;
$phiP = (32.93192) / $DEG_IN_RAD;

# Ecliptic
$T = 2;   # Centuries from 1900
$alphaP = 0;
$deltaP = ((23+27/60+8.26/3600) - 46.845/3600*$T - 0.0059/3600*$T*$T + 0.00181/3600*$T*$T*$T) / $DEG_IN_RAD;
$phiP = 0;

$alpha = 180 / $DEG_IN_RAD;
$delta = 30 / $DEG_IN_RAD;

### The transformation itself
$theta = asin( sin($delta) * cos($deltaP) - cos($delta) * sin($deltaP) * sin($alpha - $alphaP));
$y = cos($delta) * cos($deltaP) * sin($alpha - $alphaP) + sin($delta) * sin($deltaP);
$x = cos($delta) * cos($alpha - $alphaP);
$phi = atan2($y,$x) + $phiP;

$theta *= $DEG_IN_RAD;
$phi *= $DEG_IN_RAD;

print "$phi $theta\n";
