Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] Use static arrays #46

Merged
merged 1 commit into from
Oct 26, 2017
Merged

[WIP] Use static arrays #46

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions REQUIRE
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1 +1,2 @@
julia 0.6
StaticArrays
1 change: 1 addition & 0 deletions src/AstroLib.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ __precompile__()

module AstroLib
using Base.Dates
using StaticArrays

# Note on function definitions in this package. Most functions are defined as
# follows:
Expand Down
168 changes: 84 additions & 84 deletions src/baryvel.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,96 +1,96 @@
# This file is a part of AstroLib.jl. License is MIT "Expat".

const dcfel = [ 1.7400353 6.2833195099091e2 5.2796e-6 ;
6.2565836 6.2830194572674e2 -2.6180e-6 ;
4.7199666 8.3997091449254e3 -1.9780e-5 ;
0.19636505 8.4334662911720e3 -5.6044e-5 ;
4.1547339 5.2993466764997e1 5.8845e-6 ;
4.6524223 2.1354275911213e1 5.6797e-6 ;
4.2620486 7.5025342197656 5.5317e-6 ;
1.4740694 3.8377331909193 5.6093e-6 ]

const ccsel = [ 1.675104e-2 -4.179579e-5 -1.260516e-7 ;
2.220221e-1 2.809917e-2 1.852532e-5 ;
1.589963 3.418075e-2 1.430200e-5 ;
2.994089 2.590824e-2 4.155840e-6 ;
8.155457e-1 2.486352e-2 6.836840e-6 ;
1.735614 1.763719e-2 6.370440e-6 ;
1.968564 1.524020e-2 -2.517152e-6 ;
1.282417 8.703393e-3 2.289292e-5 ;
2.28082 1.918010e-2 4.484520e-6 ;
4.833473e-2 1.641773e-4 -4.654200e-7 ;
5.589232e-2 -3.455092e-4 -7.388560e-7 ;
4.634443e-2 -2.658234e-5 7.757000e-8 ;
8.997041e-3 6.329728e-6 -1.939256e-9 ;
2.284178e-2 -9.941590e-5 6.787400e-8 ;
4.350267e-2 -6.839749e-5 -2.714956e-7 ;
1.348204e-2 1.091504e-5 6.903760e-7 ;
3.106570e-2 -1.665665e-4 -1.590188e-7 ]

const dcargs = [ 5.0974222 -7.8604195454652e2 ;
3.9584962 -5.7533848094674e2 ;
1.633807 -1.1506769618935e3 ;
2.5487111 -3.9302097727326e2 ;
4.9255514 -5.8849265665348e2 ;
1.3363463 -5.5076098609303e2 ;
1.6072053 -5.2237501616674e2 ;
1.362948 -1.1790629318198e3 ;
5.5657014 -1.0977134971135e3 ;
5.0708205 -1.5774000881978e2 ;
3.9318944 5.296346478e1 ;
4.8989497 3.9809289073258e1 ;
1.3097446 7.7540959633708e1 ;
3.5147141 7.9618578146517e1 ;
3.5413158 -5.4868336758022e2 ]

const ccamps = [-2.279594e-5 1.407414e-5 8.273188e-6 1.340565e-5 -2.490817e-7 ;
-3.494537e-5 2.860401e-7 1.289448e-7 1.627237e-5 -1.823138e-7 ;
6.593466e-7 1.322572e-5 9.258695e-6 -4.674248e-7 -3.646275e-7 ;
1.140767e-5 -2.049792e-5 -4.747930e-6 -2.638763e-6 -1.245408e-7 ;
9.516893e-6 -2.748894e-6 -1.319381e-6 -4.549908e-6 -1.864821e-7 ;
7.310990e-6 -1.924710e-6 -8.772849e-7 -3.334143e-6 -1.745256e-7 ;
-2.603449e-6 7.359472e-6 3.168357e-6 1.119056e-6 -1.655307e-7 ;
-3.228859e-6 1.308997e-7 1.013137e-7 2.403899e-6 -3.736225e-7 ;
3.442177e-7 2.671323e-6 1.832858e-6 -2.394688e-7 -3.478444e-7 ;
8.702406e-6 -8.421214e-6 -1.372341e-6 -1.455234e-6 -4.998479e-8 ;
-1.488378e-6 -1.251789e-5 5.226868e-7 -2.049301e-7 0 ;
-8.043059e-6 -2.991300e-6 1.473654e-7 -3.154542e-7 0 ;
3.699128e-6 -3.316126e-6 2.901257e-7 3.407826e-7 0 ;
2.550120e-6 -1.241123e-6 9.901116e-8 2.210482e-7 0 ;
-6.351059e-7 2.341650e-6 1.061492e-6 2.878231e-7 0 ]

const ccsec = [ 1.289600e-6 5.550147e-1 2.076942 ;
3.102810e-5 4.035027 3.525565e-1 ;
9.124190e-6 9.990265e-1 2.622706 ;
9.793240e-7 5.508259 1.559103e01 ]

const dcargm = [5.1679830 8.3286911095275e3 ;
5.4913150 -7.2140632838100e3 ;
5.9598530 1.5542754389685e4 ]

const ccampm = [ 1.097594e-1 2.896773e-7 5.450474e-2 1.438491E-7 ;
-2.223581e-2 5.083103e-8 1.002548e-2 -2.291823E-8 ;
1.148966e-2 5.658888e-8 8.249439e-3 4.063015E-8 ]

const ccpamv = [8.326827e-11, 1.843484e-11, 1.988712e-12, 1.881276e-12]
const dcfel = @SMatrix [ 1.7400353 6.2833195099091e2 5.2796e-6 ;
6.2565836 6.2830194572674e2 -2.6180e-6 ;
4.7199666 8.3997091449254e3 -1.9780e-5 ;
0.19636505 8.4334662911720e3 -5.6044e-5 ;
4.1547339 5.2993466764997e1 5.8845e-6 ;
4.6524223 2.1354275911213e1 5.6797e-6 ;
4.2620486 7.5025342197656 5.5317e-6 ;
1.4740694 3.8377331909193 5.6093e-6 ]

const ccsel = @SMatrix [ 1.675104e-2 -4.179579e-5 -1.260516e-7 ;
2.220221e-1 2.809917e-2 1.852532e-5 ;
1.589963 3.418075e-2 1.430200e-5 ;
2.994089 2.590824e-2 4.155840e-6 ;
8.155457e-1 2.486352e-2 6.836840e-6 ;
1.735614 1.763719e-2 6.370440e-6 ;
1.968564 1.524020e-2 -2.517152e-6 ;
1.282417 8.703393e-3 2.289292e-5 ;
2.28082 1.918010e-2 4.484520e-6 ;
4.833473e-2 1.641773e-4 -4.654200e-7 ;
5.589232e-2 -3.455092e-4 -7.388560e-7 ;
4.634443e-2 -2.658234e-5 7.757000e-8 ;
8.997041e-3 6.329728e-6 -1.939256e-9 ;
2.284178e-2 -9.941590e-5 6.787400e-8 ;
4.350267e-2 -6.839749e-5 -2.714956e-7 ;
1.348204e-2 1.091504e-5 6.903760e-7 ;
3.106570e-2 -1.665665e-4 -1.590188e-7 ]

const dcargs = @SMatrix [ 5.0974222 -7.8604195454652e2 ;
3.9584962 -5.7533848094674e2 ;
1.633807 -1.1506769618935e3 ;
2.5487111 -3.9302097727326e2 ;
4.9255514 -5.8849265665348e2 ;
1.3363463 -5.5076098609303e2 ;
1.6072053 -5.2237501616674e2 ;
1.362948 -1.1790629318198e3 ;
5.5657014 -1.0977134971135e3 ;
5.0708205 -1.5774000881978e2 ;
3.9318944 5.296346478e1 ;
4.8989497 3.9809289073258e1 ;
1.3097446 7.7540959633708e1 ;
3.5147141 7.9618578146517e1 ;
3.5413158 -5.4868336758022e2 ]

const ccamps = @SMatrix [-2.279594e-5 1.407414e-5 8.273188e-6 1.340565e-5 -2.490817e-7 ;
-3.494537e-5 2.860401e-7 1.289448e-7 1.627237e-5 -1.823138e-7 ;
6.593466e-7 1.322572e-5 9.258695e-6 -4.674248e-7 -3.646275e-7 ;
1.140767e-5 -2.049792e-5 -4.747930e-6 -2.638763e-6 -1.245408e-7 ;
9.516893e-6 -2.748894e-6 -1.319381e-6 -4.549908e-6 -1.864821e-7 ;
7.310990e-6 -1.924710e-6 -8.772849e-7 -3.334143e-6 -1.745256e-7 ;
-2.603449e-6 7.359472e-6 3.168357e-6 1.119056e-6 -1.655307e-7 ;
-3.228859e-6 1.308997e-7 1.013137e-7 2.403899e-6 -3.736225e-7 ;
3.442177e-7 2.671323e-6 1.832858e-6 -2.394688e-7 -3.478444e-7 ;
8.702406e-6 -8.421214e-6 -1.372341e-6 -1.455234e-6 -4.998479e-8 ;
-1.488378e-6 -1.251789e-5 5.226868e-7 -2.049301e-7 0 ;
-8.043059e-6 -2.991300e-6 1.473654e-7 -3.154542e-7 0 ;
3.699128e-6 -3.316126e-6 2.901257e-7 3.407826e-7 0 ;
2.550120e-6 -1.241123e-6 9.901116e-8 2.210482e-7 0 ;
-6.351059e-7 2.341650e-6 1.061492e-6 2.878231e-7 0 ]

const ccsec = @SMatrix [ 1.289600e-6 5.550147e-1 2.076942 ;
3.102810e-5 4.035027 3.525565e-1 ;
9.124190e-6 9.990265e-1 2.622706 ;
9.793240e-7 5.508259 1.559103e01 ]

const dcargm = @SMatrix [5.1679830 8.3286911095275e3 ;
5.4913150 -7.2140632838100e3 ;
5.9598530 1.5542754389685e4 ]

const ccampm = @SMatrix [ 1.097594e-1 2.896773e-7 5.450474e-2 1.438491E-7 ;
-2.223581e-2 5.083103e-8 1.002548e-2 -2.291823E-8 ;
1.148966e-2 5.658888e-8 8.249439e-3 4.063015E-8 ]

const ccpamv = @SVector [8.326827e-11, 1.843484e-11, 1.988712e-12, 1.881276e-12]

function _baryvel(dje::T) where {T<:AbstractFloat}
# Time arguments.
dt = (dje - 2415020) / JULIANCENTURY
tvec = [1; dt; dt*dt]
tvec = SVector(1, dt, dt * dt)

temp = rem.(dcfel * tvec, 2 * T(pi))
dml, g = temp[1:2]
forbel = temp[2:8]
deps = rem.(dot(tvec, [0.4093198; -2.271110e-4; -2.860401e-8]), 2 * T(pi))
dml, g = @view temp[1:2]
forbel = @view temp[2:8]
deps = rem.(dot(tvec, SVector(0.4093198, -2.271110e-4, -2.860401e-8)), 2 * T(pi))
sorbel = rem.(ccsel * tvec, 2 * T(pi))
e = sorbel[1]

# Secular perturbations in longitude.
sn = sin.(rem.(ccsec[:,2:3] * tvec[1:2], 2 * T(pi)))
sn = sin.(rem.(@view(ccsec[:,2:3]) * @view(tvec[1:2]), 2 * T(pi)))

# Periodic perturbations of the emb (earth-moon barycenter).
pertl = dot(ccsec[:,1], sn) + dt * sn[3] * -7.757020e-8
pertl = dot(@view(ccsec[:,1]), sn) + dt * sn[3] * -7.757020e-8
pertr = zero(T)
pertld = zero(T)
pertrd = zero(T)
Expand Down Expand Up @@ -166,7 +166,7 @@ function _baryvel(dje::T) where {T<:AbstractFloat}
dyabd = dcosep * dybd - dsinep * dzbd
dzabd = dsinep * dybd + dcosep * dzbd

return [dxhd, dyahd, dzahd], [dxbd, dyabd, dzabd]
return SVector(dxhd, dyahd, dzahd), SVector(dxbd, dyabd, dzabd)
end

"""
Expand Down Expand Up @@ -216,17 +216,17 @@ baryvel(dje::Real) = _baryvel(float(dje))

function baryvel(dje::AbstractFloat)
v1, v2 = _baryvel(dje)
dvelh = AU / 1000 * v1
dvelb = AU / 1000 * v2
dvelh = AU ./ 1000 .* v1
dvelb = AU ./ 1000 .* v2
return dvelh, dvelb
end

function baryvel(dje::T, deq::T) where {T<:AbstractFloat}
v1, v2 = _baryvel(dje)
deqdat = ((dje - 2415020.313) / 365.24219572) + 1900
prema = premat(deqdat, deq, FK4 = true)
dvelh = (prema * v1) * AU / 1000
dvelb = (prema * v2) * AU / 1000
dvelh = (prema * v1) .* AU ./ 1000
dvelb = (prema * v2) .* AU ./ 1000
return dvelh, dvelb
end

Expand Down
34 changes: 17 additions & 17 deletions src/bprecess.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,18 +2,18 @@
# Copyright (C) 2016 Mosè Giordano.

const Mbprec =
reshape([+0.9999256795, -0.0111814828, -0.0048590040,
-0.000551, -0.238560, +0.435730,
+0.0111814828, +0.9999374849, -0.0000271557,
+0.238509, -0.002667, -0.008541,
+0.0048590039, -0.0000271771, +0.9999881946,
-0.435614, +0.012254, +0.002117,
-0.00000242389840, +0.00000002710544, +0.00000001177742,
+0.99990432, -0.01118145, -0.00485852,
-0.00000002710544, -0.00000242392702, +0.00000000006585,
+0.01118145, +0.99991613, -0.00002716,
-0.00000001177742, +0.00000000006585, -0.00000242404995,
+0.00485852, -0.00002717, +0.99996684], 6, 6)
SMatrix{6,6}(+0.9999256795, -0.0111814828, -0.0048590040,
-0.000551, -0.238560, +0.435730,
+0.0111814828, +0.9999374849, -0.0000271557,
+0.238509, -0.002667, -0.008541,
+0.0048590039, -0.0000271771, +0.9999881946,
-0.435614, +0.012254, +0.002117,
-0.00000242389840, +0.00000002710544, +0.00000001177742,
+0.99990432, -0.01118145, -0.00485852,
-0.00000002710544, -0.00000242392702, +0.00000000006585,
+0.01118145, +0.99991613, -0.00002716,
-0.00000001177742, +0.00000000006585, -0.00000242404995,
+0.00485852, -0.00002717, +0.99996684)

# Note: IDL version of `bprecess' changes in-place "muradec", "parallax" and
# "radvel". We don't do anything like this, but calculations are below,
Expand All @@ -23,14 +23,14 @@ function _bprecess(ra::T, dec::T, parallax::T, radvel::T,
@assert length(muradec) == 2
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
r0 = SVector(cosra*cosdec, sinra*cosdec, sindec)
r0_dot = SVector(-muradec[1]*sinra*cosdec - muradec[2]*cosra*sindec,
muradec[1]*cosra*cosdec - muradec[2]*sinra*sindec,
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_dot = @view R_1[4:6]
r1 .+= deg2rad.(r1_dot .* (epoch .- 1950) ./ 360000)
A = A_precess .+ deg2rad.(A_dot_precess .* (epoch .- 1950) ./ 360000)
else
Expand Down
4 changes: 2 additions & 2 deletions src/common.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,8 +43,8 @@ const JULIANCENTURY = 36_525 # days in one Julian century
const ct2lst_c = (280.46061837, 360.98564736629, 0.000387933, 38710000.0)

# Used in "bprecess" and "jprecess".
const A_precess = [-1.62557, -0.31919, -0.13843] ./ 1000000 # In radians
const A_dot_precess = [1.244 , -1.579, -0.660] ./ 1000 # In arc seconds per century
const A_precess = SVector(-1.62557, -0.31919, -0.13843) ./ 1000000 # In radians
const A_dot_precess = SVector(1.244 , -1.579, -0.660) ./ 1000 # In arc seconds per century

"""
List of observing sites. The observatories have `Observatory` type.
Expand Down
8 changes: 5 additions & 3 deletions src/gal_uvw.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,14 +1,16 @@
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

const A_G = SVector(0.0548755604, +0.4941094279, -0.8676661490,
0.8734370902, -0.4448296300, -0.1980763734,
0.4838350155, 0.7469822445, +0.4559837762)


function _gal_uvw(ra::T, dec::T, pmra::T, pmdec::T, vrad::T,
plx::T, lsr::Bool) where {T<:AbstractFloat}
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];
[ 0.4838350155, 0.7469822445, +0.4559837762]]
vec2 = k*pmra/plx
vec3 = k*pmdec/plx
u = ( A_G[1]*cosra*cosdec + A_G[4]*sinra*cosdec + A_G[7]*sindec)*vrad +
Expand Down
15 changes: 7 additions & 8 deletions src/geo2mag.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,18 +14,17 @@ function _geo2mag(lat::T, long::T, pole_lat::T, pole_long::T) where {T<:Abstract
# the Greenwich meridian to the meridian containing the magnetic dipole
# pole.
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]
geolong2maglong = SMatrix{3,3}(cos_pole_long, -sin_pole_long, zero(T),
sin_pole_long, cos_pole_long, zero(T),
zero(T), zero(T), one(T))

# Second rotation: in the plane of the current meridian from geographic pole
# to magnetic dipole pole.
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
tomaglat = SMatrix{3,3}( c, zero(T), s,
zero(T), one(T), zero(T),
-s, zero(T), c)
out = tomaglat * geolong2maglong * SVector(x, y, z)

maglat = rad2deg(atan2(out[3], hypot(out[1], out[2])))
maglong = rad2deg(atan2(out[2], out[1]))
Expand Down
20 changes: 10 additions & 10 deletions src/helio.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,15 +2,15 @@

# dpdt gives the time rate of change of the mean orbital quantities
# dpdt elements taken from https://ssd.jpl.nasa.gov/txt/p_elem_t1.txt
const dpd = [ 0.00000037 0.00001906 -0.00594749 -0.12534081 0.16047689 149472.67411175;
0.00000390 -0.00004107 -0.00078890 -0.27769418 0.00268329 58517.81538729;
0.00000562 -0.00004392 -0.01294668 0.0 0.32327364 35999.37244981;
0.00001847 0.00007882 -0.00813131 -0.29257343 0.44441088 19140.30268499;
-0.00011607 -0.00013253 -0.00183714 0.20469106 0.21252668 3034.74612775;
-0.00125060 -0.00050991 0.00193609 -0.28867794 -0.41897216 1222.49362201;
-0.00196176 -0.00004397 -0.00242939 0.04240589 0.40805281 428.48202785;
0.00026291 0.00005105 0.00035372 -0.00508664 -0.32241464 218.45945325;
-0.00031596 0.00005170 0.00004818 -0.01183482 -0.04062942 145.20780515]
const dpd = @SMatrix [ 0.00000037 0.00001906 -0.00594749 -0.12534081 0.16047689 149472.67411175;
0.00000390 -0.00004107 -0.00078890 -0.27769418 0.00268329 58517.81538729;
0.00000562 -0.00004392 -0.01294668 0.0 0.32327364 35999.37244981;
0.00001847 0.00007882 -0.00813131 -0.29257343 0.44441088 19140.30268499;
-0.00011607 -0.00013253 -0.00183714 0.20469106 0.21252668 3034.74612775;
-0.00125060 -0.00050991 0.00193609 -0.28867794 -0.41897216 1222.49362201;
-0.00196176 -0.00004397 -0.00242939 0.04240589 0.40805281 428.48202785;
0.00026291 0.00005105 0.00035372 -0.00508664 -0.32241464 218.45945325;
-0.00031596 0.00005170 0.00004818 -0.01183482 -0.04062942 145.20780515]

const record = Dict(1=>"mercury", 2=>"venus", 3=>"earth", 4=>"mars", 5=>"jupiter",
6=>"saturn", 7=>"uranus", 8=>"neptune", 9=>"pluto")
Expand All @@ -22,7 +22,7 @@ function _helio(jd::T, num::Integer, radians::Bool) where {T<:AbstractFloat}
end
t = (jd - J2000) / JULIANCENTURY
body = record[num]
dpdt = dpd * t
dpdt = dpd .* t
a = planets[body].axis/AU + dpdt[num, 1]
eccen = planets[body].ecc + dpdt[num, 2]
n = deg2rad(0.9856076686 / (a * sqrt(a) ))
Expand Down
Loading