# R/astro.R In tripEstimation: Metropolis Sampler and Supporting Functions for Estimating Animal Movement from Archival Tags and Satellite Fixes

````astro` <-
function(lon, lat, astro.calc) {
RA <- astro.calc\$RA
DELTA <- astro.calc\$DEC
MODJD <- astro.calc\$MJD

## beware I'm moving to cartographic longitude
TAU <- 15.0 * (LMST(MODJD, -lon) - RA)
EQUHOR(DELTA, TAU, lat)

}

`EQUHOR` <-
function(DEC,TAU,PHI) {

DEC <- DEC * (pi/180)
TAU <- TAU * (pi/180)
PHI <- PHI * (pi/180)

CS_PHI = cos(PHI); SN_PHI= sin(PHI);
CS_DEC = cos(DEC); SN_DEC =sin(DEC); CS_TAU=cos(TAU);

X =CS_DEC*SN_PHI*CS_TAU - SN_DEC*CS_PHI;
Y =CS_DEC * sin(TAU);
Z =CS_DEC*CS_PHI*CS_TAU + SN_DEC*SN_PHI;

# print(list(X, Y, Z))
HAZ <- POLAR(X,Y,Z)
#list(H = HAZ\$theta, AZ = HAZ\$phi)
HAZ
}

`FRAC` <-
function(x) {
## return the fractional portion, this is not documented in the reference
##   just relies on (assuming) how Pascal's TRUNC behaves
x <- x - trunc(x)
#if (x < 0) x <- x + 1
x
}

`julcent` <-
function(time) {
## julian centuries
tm <- as.POSIXlt(time, tz = "GMT")
hh <- tm\$hour
mm <- tm\$min
ss <- tm\$sec
jday <- julday(tm) + (hh + (mm + ss/60)/60)/24
(jday - 2451545)/36525
}

`LMST` <-
function(MJDay, LAMBDA) {
## (* LMST: local mean sidereal time
MJD0 = trunc(MJDay);

UT = (MJDay-MJD0)*24;
time =(MJD0-51544.5)/36525.0;
GMST =6.697374558 + 1.0027379093*UT +
(8640184.812866+(0.093104 - 6.2E-6*time)*time)*time/3600.0;
24.0*FRAC( (GMST-LAMBDA/15.0) / 24.0 );

}

`lunar` <-
function(time) {
ARC = 206264.8062;
P2  = 6.283185307;
COSEPS = 0.91748;
SINEPS = 0.39778;
stime <- julcent(time)
#    (* mean elements of lunar orbit *)
L0 =   FRAC(0.606433+1336.855225*stime); #(* mean longitude Moon (in rev) *)
L = P2*FRAC(0.374897+1325.552410*stime); #(* mean anomaly of the Moon     *)
LS =P2*FRAC(0.993133+  99.997361*stime); #(* mean anomaly of the Sun      *)
dd = P2*FRAC(0.827361+1236.853086*stime); #(* diff. longitude Moon-Sun     *)
ff = P2*FRAC(0.259086+1342.227825*stime); #(* mean argument of latitude    *)
DL = 22640*sin(L) - 4586*sin(L-2*dd) + 2370*sin(2*dd) +  769*sin(2*L)  -
668*sin(LS)- 412*sin(2*ff) - 212*sin(2*L-2*dd) - 206*sin(L+LS-2*dd) +
192*sin(L+2*dd) - 165*sin(LS-2*dd) - 125*sin(dd) - 110*sin(L+LS) +
148*sin(L-LS) - 55*sin(2*ff-2*dd);
S = ff + (DL+412*sin(2*ff)+541*sin(LS)) / ARC;
H = ff - 2*dd;
N = -526*sin(H) + 44*sin(L+H) - 31*sin(-L+H) - 23*sin(LS+H) +
11*sin(-LS+H) -25*sin(-2*L+ff) + 21*sin(-L+ff);
L_MOON = P2 * FRAC ( L0 + DL/1296E3 ); #(* in rad *)
B_MOON = ( 18520.0* sin(S) + N ) / ARC; #(* in rad *)
#(* equatorial coordinates *)
CB = cos(B_MOON);
X = CB*cos(L_MOON); V =CB*sin(L_MOON); W =sin(B_MOON);
Y =COSEPS*V-SINEPS*W; Z=SINEPS*V+COSEPS*W; RHO=sqrt(1.0-Z*Z);
DEC = (360.0/P2)*atan(Z/RHO);
RA  = ( 48.0/P2)*atan(Y/(X+RHO)); ifelse(RA<0, RA+24.0, RA)

MODJD <- MJD(time)
list(RA = RA, DEC = DEC, MJD = MODJD)
}

`mini.sun` <-
function(time) {

P2  = 6.283185307;
COSEPS = 0.91748;
SINEPS = 0.39778;
stime <- julcent(time)
M  = P2*FRAC(0.993133+99.997361*stime);
DL = 6893.0*sin(M)+72.0*sin(2*M);
L  = P2*FRAC(0.7859453 + M/P2 + (6191.2*stime+DL)/1296E3);
SL = sin(L);
X =cos(L);
Y = COSEPS*SL;

Z = SINEPS*SL;
RHO = sqrt(1.0-Z*Z);
DEC = (360.0/P2)*atan(Z/RHO);

RA  = ( 48.0/P2)*atan(Y/(X+RHO));
RA <- ifelse(RA < 0, RA+24.0, RA);

MODJD <- MJD(time)
list(RA = RA, DEC = DEC, MJD = MODJD)
}

`MJD` <-
function(date) {
## modified julian date (it's a smaller number)
date <- as.POSIXlt(date, tz = "GMT")
YEAR <- as.numeric(format(date, "%Y"))
MONTH <- as.numeric(format(date, "%m"))
DAY <- as.numeric(format(date, "%d"))
HOUR <- date\$hour + date\$min/60 + date\$sec/3600
A =10000.0*YEAR+100.0*MONTH+DAY;
YEAR = ifelse(MONTH <= 2, YEAR-1, YEAR)
MONTH = ifelse(MONTH <= 2, MONTH + 12, MONTH);

B = ifelse(A<=15821004.1, -2 + trunc((YEAR+4716)/4)-1179,
trunc(YEAR/400)-trunc(YEAR/100)+trunc(YEAR/4))
A =365.0*YEAR-679004.0;
A+B+trunc(30.6001*(MONTH+1))+DAY+HOUR/24.0;
}

`POLAR` <-
function(X,Y,Z) {
RHO =X*X+Y*Y;  R = sqrt(RHO+Z*Z)
PHI <- atan2(Y,X)*180/pi ; PHI <- ifelse(PHI<0, PHI+360.0, PHI);

##Thanks to Nick.Ellis@csiro.au for the report.  2009-05-29
THETA <- atan2(Z, sqrt(RHO)) * 180/pi # this line had error, was RHO now sqrt(RHO)

list(r = R, theta= THETA, phi = PHI)
}
```

## Try the tripEstimation package in your browser

Any scripts or data that you put into this service are public.

tripEstimation documentation built on May 2, 2019, 8:51 a.m.