home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
rtsi.com
/
2014.01.www.rtsi.com.tar
/
www.rtsi.com
/
OS9
/
OSK
/
EFFO
/
pd8.lzh
/
SRC
/
riset.c
< prev
next >
Wrap
C/C++ Source or Header
|
1990-04-13
|
2KB
|
85 lines
#include <stdio.h>
#include <math.h>
#include "astro.h"
/* given the true geocentric ra and dec of an object, in radians, the observer's
* latitude in radians, and a horizon displacement correction, dis, in radians
* find the local sidereal times and azimuths of rising and setting, lstr/s
* and azr/s, in radians, respectively.
* dis is the vertical displacement from the true position of the horizon. it
* is positive if the apparent position is higher than the true position.
* said another way, it is positive if the shift causes the object to spend
* longer above the horizon. for example, atmospheric refraction is typically
* assumed to produce a vertical shift of 34 arc minutes at the horizon; dis
* would then take on the value +9.89e-3 (radians). On the other hand, if
* your horizon has hills such that your apparent horizon is, say, 1 degree
* above sea level, you would allow for this by setting dis to -1.75e-2
* (radians).
* status: 0: normal; 1: never rises; -1: circumpolar; 2: trouble.
*/
riset (ra, dec, lat, dis, lstr, lsts, azr, azs, status)
double ra, dec;
double lat, dis;
double *lstr, *lsts;
double *azr, *azs;
int *status;
{
static double lastlat = 0, slat = 0, clat = 1.0;
double dec1, sdec, cdec, tdec;
double psi, spsi, cpsi;
double h, dh, ch; /* hr angle, delta and cos */
/* avoid needless sin/cos since latitude doesn't change often */
if (lat != lastlat) {
clat = cos(lat);
slat = sin(lat);
lastlat = lat;
}
/* can't cope with objects very near the celestial pole nor if we
* are located very near the earth's poles.
*/
cdec = cos(dec);
if (fabs(cdec*clat) < 1e-10) {
/* trouble */
*status = 2;
return;
}
cpsi = slat/cdec;
if (cpsi>1.0) cpsi = 1.0;
else if (cpsi<-1.0) cpsi = -1.0;
psi = acos(cpsi);
spsi = sin(psi);
dh = dis*spsi;
dec1 = dec + (dis*cpsi);
sdec = sin(dec1);
tdec = tan(dec1);
ch = (-1*slat*tdec)/clat;
if (ch < -1) {
/* circumpolar */
*status = -1;
return;
}
if (ch > 1) {
/* never rises */
*status = 1;
return;
}
*status = 0;
h = acos(ch)+dh;
*lstr = 24+radhr(ra-h);
*lsts = radhr(ra+h);
range (lstr, 24.0);
range (lsts, 24.0);
*azr = acos(sdec/clat);
range (azr, 2*PI);
*azs = 2*PI - *azr;
}