Use more sincos

JuliaAstro · Oct 12, 2017 · 3538f95 · 3538f95
1 parent 97219f5
commit 3538f95
Show file tree

Hide file tree

Showing 26 changed files with 193 additions and 167 deletions.
diff --git a/src/aitoff.jl b/src/aitoff.jl
@@ -7,9 +7,11 @@ function aitoff(l::T, b::T) where {T<:AbstractFloat}
     delta = deg2rad(b)
     r2 = sqrt(T(2))
     f = 2*r2/pi
-    cdec = cos(delta)
-    denom = sqrt(1 + cdec*cos(alpha2))
-    return rad2deg(cdec*sin(alpha2)*2*r2/denom/f), rad2deg(sin(delta)*r2/denom/f)
+    sin_alpha2, cos_alpha2 = sincos(alpha2)
+    sin_delta, cos_delta = sincos(delta)
+    denom = sqrt(1 + cos_delta * cos_alpha2)
+    return rad2deg(cos_delta * sin_alpha2 * 2 * r2 / denom / f),
+           rad2deg(sin_delta * r2 / denom / f)
 end
 
 """

diff --git a/src/altaz2hadec.jl b/src/altaz2hadec.jl
@@ -6,14 +6,17 @@ function altaz2hadec(alt::T, az::T, lat::T) where {T<:AbstractFloat}
     alt_r = deg2rad(alt)
     az_r = deg2rad(az)
     lat_r = deg2rad(lat)
+    sin_alt_r, cos_alt_r = sincos(alt_r)
+    sin_az_r, cos_az_r = sincos(az_r)
+    sin_lat_r, cos_lat_r = sincos(lat_r)
     # Find local hour angle (in degrees, from 0. to 360.).
-    ha = rad2deg(atan2(-sin(az_r)*cos(alt_r),
-                       -cos(az_r)*sin(lat_r)*cos(alt_r) +
-                       sin(alt_r)*cos(lat_r)))
+    ha = rad2deg(atan2(-sin_az_r * cos_alt_r,
+                       -cos_az_r * sin_lat_r * cos_alt_r +
+                       sin_alt_r * cos_lat_r))
     ha = mod(ha, 360)
     # Find declination (positive if north of Celestial Equator, negative if
     # south)
-    sindec = sin(lat_r)*sin(alt_r) + cos(lat_r)*cos(alt_r)*cos(az_r)
+    sindec = sin_lat_r * sin_alt_r + cos_lat_r * cos_alt_r * cos_az_r
     dec = rad2deg(asin(sindec))  # convert dec to degrees
     return ha, dec
 end

diff --git a/src/baryvel.jl b/src/baryvel.jl
@@ -151,10 +151,11 @@ function _baryvel(dje::T) where {T<:AbstractFloat}
     for i = 1:4
         plon = forbel[i+3]
         pomg = sorbel[i+1]
+        sin_pomg, cos_pomg = sincos(pomg)
         pecc = sorbel[i+9]
-        tl = rem(plon + 2 * pecc * sin(plon-pomg), 2 * T(pi))
-        dxbd += ccpamv[i] * (sin(tl) + pecc * sin(pomg))
-        dybd -= ccpamv[i] * (cos(tl) + pecc * cos(pomg))
+        sin_tl, cos_tl = sincos(rem(plon + 2 * pecc * sin(plon-pomg), 2 * T(pi)))
+        dxbd += ccpamv[i] * (sin_tl + pecc * sin_pomg)
+        dybd -= ccpamv[i] * (cos_tl + pecc * cos_pomg)
         dzbd -= ccpamv[i] * sorbel[i+13] * cos(plon - sorbel[i+5])
     end
 

diff --git a/src/bprecess.jl b/src/bprecess.jl
@@ -21,22 +21,20 @@ const Mbprec =
 function _bprecess(ra::T, dec::T, parallax::T, radvel::T,
                    epoch::T, muradec::Vector{T}) where {T<:AbstractFloat}
     @assert length(muradec) == 2
-    cosra  = cosd(ra)
-    sinra  = sind(ra)
-    cosdec = cosd(dec)
-    sindec = sind(dec)
-    ra1950 = dec1950 = zero(T)
-    A  = copy(A_precess)
+    sinra,  cosra  = sincos(deg2rad(ra))
+    sindec, cosdec = sincos(deg2rad(dec))
     r0 = [cosra*cosdec,  sinra*cosdec,  sindec]
     r0_dot = [-muradec[1]*sinra*cosdec - muradec[2]*cosra*sindec,
                muradec[1]*cosra*cosdec - muradec[2]*sinra*sindec,
-               muradec[2]*cosdec] + 21.095*radvel*parallax*r0
+               muradec[2]*cosdec] .+ 21.095 .* radvel * parallax * r0
     R_1 = Mbprec * vcat(r0, r0_dot)
     r1 = R_1[1:3]
     r1_dot = R_1[4:6]
     if isfinite(epoch)
         r1 .+= deg2rad.(r1_dot .* (epoch .- 1950) ./ 360000)
-        A  .+= deg2rad.(A_dot_precess .* (epoch .- 1950) ./ 360000)
+        A = A_precess .+ deg2rad.(A_dot_precess .* (epoch .- 1950) ./ 360000)
+    else
+        A = A_precess
     end
     rmag = vecnorm(r1)
     s1 = r1 ./ rmag

diff --git a/src/co_nutate.jl b/src/co_nutate.jl
@@ -4,9 +4,11 @@ function _co_nutate(jd::T, ra::T, dec::T) where {T<:AbstractFloat}
     d_psi, d_eps = nutate(jd)
     eps = mean_obliquity(jd) + sec2rad(d_eps)
     se, ce = sincos(eps)
-    x = cosd(ra) * cosd(dec)
-    y = sind(ra) * cosd(dec)
-    z = sind(dec)
+    sin_ra, cos_ra = sincos(deg2rad(ra))
+    sin_dec, cos_dec = sincos(deg2rad(dec))
+    x = cos_ra * cos_dec
+    y = sin_ra * cos_dec
+    z = sin_dec
     x2 = x - sec2rad(y * ce + z * se) * d_psi
     y2 = y + sec2rad(x * ce * d_psi - z * d_eps)
     z2 = z + sec2rad(x * se * d_psi + y * d_eps)

diff --git a/src/euler.jl b/src/euler.jl
@@ -22,15 +22,13 @@ function _euler(ai::T, bi::T, select::Integer, FK4::Bool, radians::Bool) where {
         phi = (4.9368292465, 0.574770433, 0.0, 0.0, 4.71279419371, 0.11142137093)
     end
 
-    if radians
-        ao = ai - phi[select]
-        cb = cos(bi)
-        x = (cb*sin(ao), cb*cos(ao), sin(bi))
-    else
-        ao = deg2rad(ai) - phi[select]
-        cb = cosd(bi)
-        x = (cb*sin(ao), cb*cos(ao), sind(bi))
+    if !radians
+        ai = deg2rad(ai)
+        bi = deg2rad(bi)
     end
+    sa, ca = sincos(ai - phi[select])
+    sb, cb = sincos(bi)
+    x = (cb * sa, cb * ca, sb)
     bo = ctheta[select]*x[3] - stheta[select]*x[1]
     ao = mod2pi(atan2(ctheta[select]*x[1] + stheta[select]*x[3], x[2]) + psi[select])
     bo = asin(bo)

diff --git a/src/gal_uvw.jl b/src/gal_uvw.jl
@@ -3,10 +3,8 @@
 
 function _gal_uvw(ra::T, dec::T, pmra::T, pmdec::T, vrad::T,
                   plx::T, lsr::Bool) where {T<:AbstractFloat}
-    cosdec = cosd(dec)
-    sindec = sind(dec)
-    cosra  = cosd(ra)
-    sinra  = sind(ra)
+    sindec, cosdec = sincos(deg2rad(dec))
+    sinra, cosra = sincos(deg2rad(ra))
     k = 4.740470463533348 # = 149597870.7/(86400*365.25) = 1 AU/year in km/s
     A_G = [[ 0.0548755604, +0.4941094279, -0.8676661490];
            [ 0.8734370902, -0.4448296300, -0.1980763734];

diff --git a/src/geo2eci.jl b/src/geo2eci.jl
@@ -3,14 +3,15 @@
 
 function geo2eci(lat::T, long::T, alt::T, jd::T) where {T<:AbstractFloat}
     Re    = planets["earth"].eqradius / 1000
-    lat   = deg2rad(lat)
+    sin_lat, cos_lat = sincos(deg2rad(lat))
     long  = deg2rad(long)
     gst   = ct2lst(zero(T), jd)
     sid_angle = gst*pi/12 # Sidereal angle.
     theta = long + sid_angle # Azimuth
     altRe = alt + Re
-    r     = altRe*cos(lat)
-    return r*cos(theta), r*sin(theta), altRe*sin(lat)
+    r     = altRe * cos_lat
+    sin_theta, cos_theta = sincos(theta)
+    return r * cos_theta, r * sin_theta, altRe * sin_lat
 end
 
 """

diff --git a/src/geo2geodetic.jl b/src/geo2geodetic.jl
@@ -3,11 +3,11 @@
 
 function geo2geodetic(lat::T, long::T, alt::T, eqrad::T, polrad::T) where {T<:AbstractFloat}
     e = sqrt(eqrad^2 - polrad^2)/eqrad
-    lat = deg2rad(lat)
-    long_rad = deg2rad(long)
-    x = (eqrad + alt)*cos(lat)*cos(long_rad)
-    y = (eqrad + alt)*cos(lat)*sin(long_rad)
-    z = (eqrad + alt)*sin(lat)
+    sin_lat, cos_lat = sincos(deg2rad(lat))
+    sin_long_rad, cos_long_rad = sincos(deg2rad(long))
+    x = (eqrad + alt) * cos_lat * cos_long_rad
+    y = (eqrad + alt) * cos_lat * sin_long_rad
+    z = (eqrad + alt) * sin_lat
     r = hypot(x, y)
     s    = hypot(r, z)*(1 - eqrad*sqrt((1 - e^2)/((1 - e^2)*r^2 + z^2)))
     t0   = 1 + s*sqrt(1 - (e*z)^2/(r^2 + z^2))/eqrad

diff --git a/src/geo2mag.jl b/src/geo2mag.jl
@@ -4,27 +4,27 @@
 function _geo2mag(lat::T, long::T, pole_lat::T, pole_long::T) where {T<:AbstractFloat}
     r    = 1 # Distance from planet center.  Value unimportant -- just need a length for
              # conversion to rectangular coordinates
-    lat  = deg2rad(lat)
-    long = deg2rad(long)
-    x    = r * cos(lat) * cos(long)
-    y    = r * cos(lat) * sin(long)
-    z    = r * sin(lat)
+    sin_lat, cos_lat = sincos(deg2rad(lat))
+    sin_long, cos_long = sincos(deg2rad(long))
+    x    = r * cos_lat * cos_long
+    y    = r * cos_lat * sin_long
+    z    = r * sin_lat
 
     # Compute first rotation matrix: rotation around plane of the equator, from
     # the Greenwich meridian to the meridian containing the magnetic dipole
     # pole.
-    geolong2maglong = Array{T}(3, 3)
-    geolong2maglong[:,1] = [cos(pole_long), -sin(pole_long), 0.0]
-    geolong2maglong[:,2] = [sin(pole_long),  cos(pole_long), 0.0]
-    geolong2maglong[:,3] = [           0.0,             0.0, 1.0]
+    sin_pole_long, cos_pole_long = sincos(pole_long)
+    geolong2maglong = [ cos_pole_long sin_pole_long zero(T);
+                       -sin_pole_long cos_pole_long zero(T);
+                        zero(T)       zero(T)       one(T)]
     out = geolong2maglong * [x, y, z]
 
     # Second rotation: in the plane of the current meridian from geographic pole
     # to magnetic dipole pole.
-    tomaglat = Array{T}(3, 3)
-    tomaglat[:,1] = [ cos(pi/2 - pole_lat), 0.0, sin(pi/2 - pole_lat)]
-    tomaglat[:,2] = [                  0.0, 1.0,                  0.0]
-    tomaglat[:,3] = [-sin(pi/2 - pole_lat), 0.0, cos(pi/2 - pole_lat)]
+    s, c = sincos(T(pi) / 2 - pole_lat)
+    tomaglat = [c       zero(T) -s;
+                zero(T) one(T)  zero(T);
+                s       zero(T) c]
     out = tomaglat * out
 
     maglat  = rad2deg(atan2(out[3], hypot(out[1], out[2])))

diff --git a/src/geodetic2geo.jl b/src/geodetic2geo.jl
@@ -3,10 +3,10 @@
 
 function geodetic2geo(lat::T, long::T, alt::T, eqrad::T, polrad::T) where {T<:AbstractFloat}
     e = sqrt(eqrad^2 - polrad^2)/eqrad
-    lat = deg2rad(lat)
-    beta = sqrt(1 - (e*sin(lat))^2)
-    r = (eqrad/beta + alt)*cos(lat)
-    z = (eqrad*(1 - e^2)/beta + alt)*sin(lat)
+    sin_lat, cos_lat = sincos(deg2rad(lat))
+    beta = sqrt(1 - (e * sin_lat) ^ 2)
+    r = (eqrad / beta + alt) * cos_lat
+    z = (eqrad*(1 - e ^ 2) / beta + alt) * sin_lat
     return rad2deg(atan2(z,r)), long, hypot(r, z) - eqrad
 end
 

diff --git a/src/hadec2altaz.jl b/src/hadec2altaz.jl
@@ -2,12 +2,9 @@
 # Copyright (C) 2016 Mosè Giordano.
 
 function _hadec2altaz(ha::T, dec::T, lat::T, ws::Bool) where {T<:AbstractFloat}
-    sh = sind(ha)
-    ch = cosd(ha)
-    sd = sind(dec)
-    cd = cosd(dec)
-    sl = sind(lat)
-    cl = cosd(lat)
+    sh, ch = sincos(deg2rad(ha))
+    sd, cd = sincos(deg2rad(dec))
+    sl, cl = sincos(deg2rad(lat))
 
     x = -ch*cd*sl + sd*cl
     y = -sh*cd

diff --git a/src/helio_jd.jl b/src/helio_jd.jl
@@ -10,11 +10,11 @@ function _helio_jd(date::T, ra::T, dec::T, B1950::Bool, diff::Bool) where {T<:Ab
     delta_t = (date - 33282.42345905) / JULIANCENTURY
     epsilon_sec = @evalpoly(delta_t, 44.836, -46.8495, -0.00429, 0.00181)
     epsilon = deg2rad(23.433333 + epsilon_sec/3600)
-    ra = deg2rad(ra)
-    dec = deg2rad(dec)
+    sin_ra, cos_ra = sincos(deg2rad(ra))
+    sin_dec, cos_dec = sincos(deg2rad(dec))
     x, y, z = xyz(date)
-    time = -499.00522*(cos(dec)*cos(ra)*x +
-                       (tan(epsilon)*sin(dec) + cos(dec)*sin(ra))*y)
+    time = -499.00522*(cos_dec * cos_ra * x +
+                       (tan(epsilon) * sin_dec + cos_dec * sin_ra) * y)
     if diff
         return time
     else

diff --git a/src/jprecess.jl b/src/jprecess.jl
@@ -21,14 +21,12 @@ const Mjprec =
 function _jprecess(ra::T, dec::T, parallax::T, radvel::T, epoch::T,
                    muradec::Vector{T}) where {T<:AbstractFloat}
     @assert length(muradec) == 2
-    cosra  = cosd(ra)
-    sinra  = sind(ra)
-    cosdec = cosd(dec)
-    sindec = sind(dec)
-    ra2000 = dec2000 = zero(T)
-    A  = copy(A_precess)
+    sinra, cosra  = sincos(deg2rad(ra))
+    sindec, cosdec = sincos(deg2rad(dec))
     if isfinite(epoch) && epoch != 1950
-        A  .+= deg2rad.(A_dot_precess .* (epoch .- 1950) ./ 360000)
+        A = A_precess .+ deg2rad.(A_dot_precess .* (epoch .- 1950) ./ 360000)
+    else
+        A = A_precess
     end
     r0 = [cosra*cosdec,  sinra*cosdec,  sindec]
     r0_dot = [-muradec[1]*sinra*cosdec - muradec[2]*cosra*sindec,

diff --git a/src/mag2geo.jl b/src/mag2geo.jl
@@ -4,32 +4,32 @@
 function _mag2geo(lat::T, long::T, pole_lat::T, pole_long::T) where {T<:AbstractFloat}
     r    = 1 # Distance from planet center.  Value unimportant -- just need
              # a length for conversion to rectangular coordinates
-    lat  = deg2rad(lat)
-    long = deg2rad(long)
+    sin_lat, cos_lat = sincos(deg2rad(lat))
+    sin_long, cos_long = sincos(deg2rad(long))
 
     # convert to rectangular coordinates
     #   x-axis: defined by the vector going from Earth's center towards
     #        the intersection of the equator and Greenwich's meridian.
     #   z-axis: axis of the geographic poles
     #   y-axis: defined by y=z^x
-    x = r * cos(lat) * cos(long)
-    y = r * cos(lat) * sin(long)
-    z = r * sin(lat)
+    x = r * cos_lat * cos_long
+    y = r * cos_lat * sin_long
+    z = r * sin_lat
 
     # First rotation: in the plane of the current meridian from magnetic pole to
     # geographic pole.
-    togeolat = Array{T}(3, 3)
-    togeolat[:,1] = [cos(pi/2 - pole_lat), 0.0, -sin(pi/2 - pole_lat)]
-    togeolat[:,2] = [                 0.0, 1.0,                   0.0]
-    togeolat[:,3] = [sin(pi/2 - pole_lat), 0.0,  cos(pi/2 - pole_lat)]
+    s, c = sincos(T(pi) / 2 - pole_lat)
+    togeolat = [c       zero(T) s;
+                zero(T) one(T)  zero(T);
+                -s      zero(T) c]
     out = togeolat * [x, y, z]
 
     # Second rotation matrix: rotation around plane of the equator, from the
     # meridian containing the magnetic poles to the Greenwich meridian.
-    maglong2geolong = Array{T}(3, 3)
-    maglong2geolong[:,1] = [ cos(pole_long), sin(pole_long), 0.0]
-    maglong2geolong[:,2] = [-sin(pole_long), cos(pole_long), 0.0]
-    maglong2geolong[:,3] = [            0.0,            0.0, 1.0]
+    sin_pole_long, cos_pole_long = sincos(pole_long)
+    maglong2geolong = [cos_pole_long -sin_pole_long zero(T);
+                       sin_pole_long  cos_pole_long zero(T);
+                       zero(T)        zero(T)       one(T)]
     out = maglong2geolong * out
 
     geolat  = rad2deg(atan2(out[3], hypot(out[1], out[2])))

diff --git a/src/moonpos.jl b/src/moonpos.jl
@@ -117,8 +117,11 @@ function _moonpos(jd::AbstractFloat, radians::Bool)
     ɛ = ten(23, 26) + @evalpoly(t/100, 21.448, -4680.93, -1.55, 1999.25, -51.38,
                                 -249.67, -39.05, 7.12, 27.87, 5.79, 2.45)/3600
     ɛ = deg2rad(ɛ + elong/3600)
-    ra = mod2pi(atan2(sin(λ)*cos(ɛ) - tan(β)*sin(ɛ), cos(λ)))
-    dec = asin(sin(β)*cos(ɛ) + cos(β)*sin(ɛ)*sin(λ))
+    sin_β, cos_β = sincos(β)
+    sin_ɛ, cos_ɛ = sincos(ɛ)
+    sin_λ, cos_λ = sincos(λ)
+    ra = mod2pi(atan2(sin_λ * cos_ɛ - tan(β) * sin_ɛ, cos_λ))
+    dec = asin(sin_β * cos_ɛ + cos_β * sin_ɛ * sin_λ)
     if radians
         return ra, dec, dis, λ, β
     else

diff --git a/src/mphase.jl b/src/mphase.jl
@@ -6,9 +6,12 @@ function mphase(jd::AbstractFloat)
     ras, decs = sunpos(jd, radians=true)
     # phi: geocentric elongation of the Moon from the Sun
     # inc: selenocentric (Moon centered) elongation of the Earth from the Sun
-    phi = acos(sin(decs)*sin(decm) + cos(decs)*cos(decm)*cos(ras - ram))
+    sin_decs, cos_decs = sincos(decs)
+    sin_decm, cos_decm = sincos(decm)
+    phi = acos(sin_decs * sin_decm + cos_decs * cos_decm * cos(ras - ram))
     # "dism" is in kilometers, AU in meters
-    inc = atan2(AU*sin(phi), dism*1e3 - AU*cos(phi))
+    sin_phi, cos_phi = sincos(phi)
+    inc = atan2(AU * sin_phi, dism * 1000 - AU * cos_phi)
     return (1 + cos(inc))/2
 end
 

diff --git a/src/planet_coords.jl b/src/planet_coords.jl
@@ -4,9 +4,13 @@
 function _planet_coords(date::T, num::Integer) where {T<:AbstractFloat}
     rad, long, lat = helio(date, num, true)
     rade, longe, late = helio(date, 3, true)
-    x =  rad * cos(lat) * cos(long) - rade * cos(late) * cos(longe)
-    y = rad * cos(lat) * sin(long) - rade * cos(late) * sin(longe)
-    z = rad * sin(lat) - rade * sin(late)
+    sin_lat, cos_lat = sincos(lat)
+    sin_long, cos_long = sincos(long)
+    sin_late, cos_late = sincos(late)
+    sin_longe, cos_longe = sincos(longe)
+    x = rad * cos_lat * cos_long - rade * cos_late * cos_longe
+    y = rad * cos_lat * sin_long - rade * cos_late * sin_longe
+    z = rad * sin_lat - rade * sin_late
     lamda = rad2deg(atan2(y,x))
     beta = rad2deg(atan2(z, hypot(x,y)))
     ra, dec = euler(lamda, beta, 4)

diff --git a/src/polrec.jl b/src/polrec.jl
@@ -3,10 +3,10 @@
 
 function _polrec(radius::T, angle::T, degrees::Bool) where {T<:AbstractFloat}
     if degrees
-        return radius*cos(deg2rad(angle)), radius*sin(deg2rad(angle))
-    else
-        return radius*cos(angle), radius*sin(angle)
+        angle = deg2rad(angle)
     end
+    s, c = sincos(angle)
+    return radius * c, radius * s
 end
 
 """

diff --git a/src/posang.jl b/src/posang.jl
@@ -25,9 +25,10 @@ function posang(units::Integer, ra1::T, dec1::T, ra2::T, dec2::T) where {T<:Abst
         # In any other case throw an error.
         error("units must be 0 (radians), 1 (hours, degrees) or 2 (degrees)")
     end
-    radif = ra2_rad - ra1_rad
-    angle = atan2(sin(radif), cos(dec1_rad)*tan(dec2_rad) -
-                  sin(dec1_rad)*cos(radif))
+    sin_radif, cos_radif = sincos(ra2_rad - ra1_rad)
+    sin_dec1, cos_dec1 = sincos(dec1_rad)
+    angle = atan2(sin_radif, cos_dec1 * tan(dec2_rad) -
+                  sin_dec1 * cos_radif)
     if units == 0
         return angle
     else