source: trunk/MagicSoft/slalib/mapqk.c@ 18162

Last change on this file since 18162 was 731, checked in by tbretz, 24 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 4.3 KB
Line 
1#include "slalib.h"
2#include "slamac.h"
3void slaMapqk ( double rm, double dm, double pr, double pd,
4 double px, double rv, double amprms[21],
5 double *ra, double *da )
6/*
7** - - - - - - - - -
8** s l a M a p q k
9** - - - - - - - - -
10**
11** Quick mean to apparent place: transform a star RA,Dec from
12** mean place to geocentric apparent place, given the
13** star-independent parameters.
14**
15** Use of this routine is appropriate when efficiency is important
16** and where many star positions, all referred to the same equator
17** and equinox, are to be transformed for one epoch. The
18** star-independent parameters can be obtained by calling the
19** slaMappa routine.
20**
21** If the parallax and proper motions are zero the slaMapqkz
22** routine can be used instead.
23**
24** The reference frames and timescales used are post IAU 1976.
25**
26** Given:
27** rm,dm double mean RA,Dec (rad)
28** pr,pd double proper motions: RA,Dec changes per Julian year
29** px double parallax (arcsec)
30** rv double radial velocity (km/sec, +ve if receding)
31**
32** amprms double[21] star-independent mean-to-apparent parameters:
33**
34** (0) time interval for proper motion (Julian years)
35** (1-3) barycentric position of the Earth (AU)
36** (4-6) heliocentric direction of the Earth (unit vector)
37** (7) (grav rad Sun)*2/(Sun-Earth distance)
38** (8-10) barycentric Earth velocity in units of c
39** (11) sqrt(1-v**2) where v=modulus(abv)
40** (12-20) precession/nutation (3,3) matrix
41**
42** Returned:
43** *ra,*da double apparent RA,Dec (rad)
44**
45** References:
46** 1984 Astronomical Almanac, pp B39-B41.
47** (also Lederle & Schwan, Astron. Astrophys. 134, 1-6, 1984)
48**
49** Notes:
50**
51** 1) The vectors amprms(1-3) and amprms(4-6) are referred to
52** the mean equinox and equator of epoch eq.
53**
54** 2) Strictly speaking, the routine is not valid for solar-system
55** sources, though the error will usually be extremely small.
56** However, to prevent gross errors in the case where the
57** position of the Sun is specified, the gravitational
58** deflection term is restrained within about 920 arcsec of the
59** centre of the Sun's disc. The term has a maximum value of
60** about 1.85 arcsec at this radius, and decreases to zero as
61** the centre of the disc is approached.
62**
63** Called:
64** slaDcs2c spherical to Cartesian
65** slaDvdv dot product
66** slaDmxv matrix x vector
67** slaDcc2s Cartesian to spherical
68** slaDranrm normalize angle 0-2pi
69**
70** Defined in slamac.h: DAS2R
71**
72** Last revision: 15 January 2000
73**
74** Copyright P.T.Wallace. All rights reserved.
75*/
76
77#define VF 0.21094502 /* Km/s to AU/year */
78
79{
80 int i;
81 double pmt, gr2e, ab1, eb[3], ehn[3], abv[3],
82 q[3], pxr, w, em[3], p[3], pn[3], pde, pdep1,
83 p1[3], p1dv, p2[3], p3[3];
84
85/* Unpack scalar and vector parameters */
86 pmt = amprms[0];
87 gr2e = amprms[7];
88 ab1 = amprms[11];
89 for ( i = 0; i < 3; i++ )
90 {
91 eb[i] = amprms[i+1];
92 ehn[i] = amprms[i+4];
93 abv[i] = amprms[i+8];
94 }
95
96/* Spherical to x,y,z */
97 slaDcs2c ( rm, dm, q );
98
99/* Space motion (radians per year) */
100 pxr = px * DAS2R;
101 w = VF * rv * pxr;
102 em[0] = (-pr * q[1]) - ( pd * cos ( rm ) * sin ( dm ) ) + ( w * q[0] );
103 em[1] = ( pr * q[0]) - ( pd * sin ( rm ) * sin ( dm ) ) + ( w * q[1] );
104 em[2] = ( pd * cos ( dm ) ) + ( w * q[2] );
105
106/* Geocentric direction of star (normalized) */
107 for ( i = 0; i < 3; i++ ) {
108 p[i] = q[i] + ( pmt * em[i] ) - ( pxr * eb[i] );
109 }
110 slaDvn ( p, pn, &w );
111
112/* Light deflection (restrained within the Sun's disc) */
113 pde = slaDvdv ( pn, ehn );
114 pdep1 = 1.0 + pde;
115 w = gr2e / gmax ( pdep1, 1.0e-5 );
116 for ( i = 0; i < 3; i++ ) {
117 p1[i] = pn[i] + ( w * ( ehn[i] - pde * pn[i] ) );
118 }
119
120/* Aberration (normalization omitted) */
121 p1dv = slaDvdv ( p1, abv );
122 w = 1.0 + p1dv / ( ab1 + 1.0 );
123 for ( i = 0; i < 3; i++ ) {
124 p2[i] = ab1 * p1[i] + w * abv[i];
125 }
126
127/* Precession and nutation */
128 slaDmxv ( (double(*)[3]) &amprms[12], p2, p3 );
129
130/* Geocentric apparent RA,dec */
131 slaDcc2s ( p3, ra, da );
132
133 *ra = slaDranrm ( *ra );
134}
Note: See TracBrowser for help on using the repository browser.