source: trunk/MagicSoft/slalib/pda2h.c@ 733

Last change on this file since 733 was 731, checked in by tbretz, 24 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 2.3 KB
Line 
1#include "slalib.h"
2#include "slamac.h"
3void slaPda2h ( double p, double d, double a,
4 double *h1, int *j1, double *h2, int *j2 )
5/*
6** - - - - - - - - -
7** s l a P d a 2 h
8** - - - - - - - - -
9**
10** Hour Angle corresponding to a given azimuth
11**
12** (double precision)
13**
14** Given:
15** p double latitude
16** d double declination
17** a double azimuth
18**
19** Returned:
20** *h1 double hour angle: first solution if any
21** *j1 int flag: 0 = solution 1 is valid
22** *h2 double hour angle: first solution if any
23** *j2 int flag: 0 = solution 2 is valid
24**
25** Defined in slamac.h: DPI, DPIBY2
26**
27** Last revision: 24 November 1994
28**
29** Copyright P.T.Wallace. All rights reserved.
30*/
31
32#define TINY 1e-12 /* Zone of avoidance around critical angles */
33
34{
35 double pn, an, dn, sa, ca, sasp, qt, qb, hpt, t;
36
37/* Preset status flags to OK */
38 *j1 = 0;
39 *j2 = 0;
40
41/* Adjust latitude, azimuth, declination to avoid critical values */
42 pn = slaDrange ( p );
43 if ( fabs ( fabs ( pn ) - DPIBY2 ) < TINY ) {
44 pn -= dsign ( TINY, pn);
45 } else if ( fabs ( pn ) < TINY ) {
46 pn = TINY;
47 }
48 an = slaDrange ( a );
49 if ( fabs ( fabs ( an ) - DPI ) < TINY ) {
50 an -= dsign ( TINY, an );
51 } else if ( fabs ( an ) < TINY ) {
52 an = TINY;
53 }
54 dn = slaDrange ( d );
55 if ( fabs ( fabs ( dn ) - fabs ( p ) ) < TINY ) {
56 dn -= dsign ( TINY, dn );
57 } else if ( fabs ( fabs ( dn ) - DPIBY2 ) < TINY ) {
58 dn -= dsign ( TINY, dn );
59 } else if ( fabs ( dn ) < TINY ) {
60 dn = TINY;
61 }
62
63/* Useful functions */
64 sa = sin ( an );
65 ca = cos ( an );
66 sasp = sa * sin ( pn );
67
68/* Quotient giving sin(h+t) */
69 qt = sin ( dn ) * sa * cos ( pn );
70 qb = cos ( dn ) * sqrt ( ca * ca + sasp * sasp );
71
72/* Any solutions? */
73 if ( fabs ( qt ) <= qb ) {
74
75 /* Yes: find h+t and t */
76 hpt = asin ( qt / qb );
77 t = atan2 ( sasp, - ca );
78
79 /* The two solutions */
80 *h1 = slaDrange ( hpt - t );
81 *h2 = slaDrange ( - hpt - ( t + DPI ) );
82
83 /* Reject unless h and A different signs */
84 if ( *h1 * an > 0.0 ) *j1 = - 1;
85 if ( *h2 * an > 0.0 ) *j2 = - 1;
86 } else {
87 *j1 = - 1;
88 *j2 = - 1;
89 }
90}
Note: See TracBrowser for help on using the repository browser.