aboutsummaryrefslogtreecommitdiff
path: root/src/slalib/unpcd.f
blob: 95026712e3edcbb720e43eb9a7ca360e5bc75691 (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
SUBROUTINE sla_UNPCD (DISCO,X,Y)
*+
*     - - - - - -
*      U N P C D
*     - - - - - -
*
*  Remove pincushion/barrel distortion from a distorted [x,y]
*  to give tangent-plane [x,y].
*
*  Given:
*     DISCO    d      pincushion/barrel distortion coefficient
*     X,Y      d      distorted coordinates
*
*  Returned:
*     X,Y      d      tangent-plane coordinates
*
*  Notes:
*
*  1)  The distortion is of the form RP = R*(1 + C*R**2), where R is
*      the radial distance from the tangent point, C is the DISCO
*      argument, and RP is the radial distance in the presence of
*      the distortion.
*
*  2)  For pincushion distortion, C is +ve;  for barrel distortion,
*      C is -ve.
*
*  3)  For X,Y in "radians" - units of one projection radius,
*      which in the case of a photograph is the focal length of
*      the camera - the following DISCO values apply:
*
*          Geometry          DISCO
*
*          astrograph         0.0
*          Schmidt           -0.3333
*          AAT PF doublet  +147.069
*          AAT PF triplet  +178.585
*          AAT f/8          +21.20
*          JKT f/8          +13.32
*
*  4)  The present routine is an approximate inverse to the
*      companion routine sla_PCD, obtained from two iterations
*      of Newton's method.  The mismatch between the sla_PCD and
*      sla_UNPCD routines is negligible for astrometric applications;
*      to reach 1 milliarcsec at the edge of the AAT triplet or
*      Schmidt field would require field diameters of 2.4 degrees
*      and 42 degrees respectively.
*
*  P.T.Wallace   Starlink   1 August 1994
*
*  Copyright (C) 1995 Rutherford Appleton Laboratory
*-

      IMPLICIT NONE

      DOUBLE PRECISION DISCO,X,Y

      DOUBLE PRECISION CR2,A,CR2A2,F



      CR2=DISCO*(X*X+Y*Y)
      A=(2D0*CR2+1D0)/(3D0*CR2+1D0)
      CR2A2=CR2*A*A
      F=(2D0*CR2A2*A+1D0)/(3D0*CR2A2+1D0)
      X=X*F
      Y=Y*F

      END