IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Opened 22 years ago

Closed 21 years ago

Last modified 21 years ago

#260 closed defect (fixed)

psDeproject()

Reported by: gusciora@… Owned by: Paul Price
Priority: high Milestone:
Component: PSLib ADD Version: unspecified
Severity: normal Keywords:
Cc:

Description

I'm hoping you can provide an algorithm, or formula. for the Gnomonic and
Orthographic deprojections. My attempts at deriving them from the information
in the ADD have failed: if I project a spherical coordinate onto a plane, then
deproject that planar coordinate, I don't receive the original spherical
coordinate. Also, it will be helpful if you point out which greek symbols
correspond to which spherical coordinate members of the psSphere structure;
I think I know, but I'm not 100% sure.

Attachments (1)

projections.pl (1.0 KB ) - added by Paul Price 21 years ago.
Perl code that projects and deprojects

Download all attachments as: .zip

Change History (5)

comment:1 by gusciora@…, 22 years ago

The current status is this: I couldn't get the ADD equations to work, but it's
quite possible that I was putting them together in the wrong way. I wanted
something to work for the current release and since you folks were in San Diego
last week, I looked on the web for gnomonic projections. I implemented the
formulas that are at this page:

http://mathworld.wolfram.com/GnomonicProjection.html

They basically work for several cases, but don't work for others. I'm thinking
it might be a problem with angles outside the first quadrant and the atan()
function. That web page specifically mentioned that the 2-argument arctan
function should be used, but I could not figure out how integrate it.

Personally, I don't have any preference to using the mathworks formulas over the
ADD formulas: I just wanted something to work for this release since a few other
functions are dependent on them. I will need help with either getting the ADD
equations to work, or the mathworks equations.

comment:2 by Paul Price, 21 years ago

Status: newassigned

Are you sure? It seems to work for me:

alpha,delta: 0.785398163397448 0.785398163397448
phi,theta: -0.615479708670387 0.523598775598299
x,y: -1 -1.41421356237309
phi,theta: 2.52611294491941 0.523598775598299
alpha,delta: 0.785398163397448 0.785398163397448

I will attach the code.

by Paul Price, 21 years ago

Attachment: projections.pl added

Perl code that projects and deprojects

comment:3 by Eric.VanAlst@…, 21 years ago

Code has been implemented to perform the TAN, SIN, AIT and PAR projections
specified in the ADD. The equations for the reverse of the AIT did not produce
the original. The document "Representations of celestial coordinates in FITS"
by Calabretta & Greisen shows the reverse different which has been implemented
and provide the expected results. Which equations should we use for the reverse
AIT?

ADD:

phi = 2 * arctan( 2*z*z - 1, xz)
z = sqrt( (1 - (x/2)*(x/2) - y*y) )

Calbretta & Greisen:

phi = 2* arctan( (2 * ( 2*z*z - 1)), xz)
z = sqrt ( 1 - (x/4)*(x/4) - (y/2)*(y/2) )

comment:4 by Paul Price, 21 years ago

Resolution: fixed
Status: assignedclosed

I think the problem is in getting the correct values for alphaP,deltaP. deltaP
needs to be 90 degrees so that theta has the correct range. Not sure why alphaP
needs to have 180 degrees subtracted, but it works....

The following Perl code demonstrates:

# Aitoff (AIT) projection
$deltaP = 90*pi/180;
$alphaP -= 180*pi/180;
$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));
$squiggle = sqrt(2/(1 + cos($theta)*cos($phi/2)));
$x = 2.0 * $squiggle * cos($theta) * sin($phi/2);
$y = $squiggle * sin($theta);
# Inverse Aitoff projection
$z = sqrt(1 - ($x/4)2 - ($y/2)2);
$phiNew = 2.0*atan2(2*$z2 - 1, $x*$z/2);
$thetaNew = asin($y*$z);
$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));
if ($alphaNew < 0) { $alphaNew += 2*pi; }

print "alpha,delta: $alpha $delta\n";
print "phi,theta: $phi $theta\n";
print "x,y: $x $y\n";
print "phi,theta: $phiNew $thetaNew\n";
print "alpha,delta: $alphaNew $deltaNew\n";

===>

alpha,delta: 4.71238898038469 -0.523598775598299
phi,theta: -1.22460635382238e-16 -0.523598775598299
x,y: -1.0979520198983e-16 -0.517638090205041
phi,theta: 3.14159265358979 -0.523598775598299
alpha,delta: 4.71238898038469 -0.523598775598299

Note: See TracTickets for help on using tickets.