diff --git a/src/adstring.jl b/src/adstring.jl
index 867cd85..e0098b3 100644
--- a/src/adstring.jl
+++ b/src/adstring.jl
@@ -72,7 +72,7 @@ julia> adstring(30.4, -1.23, truncate=true)
 " 02 01 35.9  -01 13 48"
 
 julia> adstring.([30.4, -15.63], [-1.23, 48.41], precision=1)
-2-element Array{String,1}:
+2-element Vector{String}:
  " 02 01 36.00  -01 13 48.0"
  " 22 57 28.80  +48 24 36.0"
 ```
diff --git a/src/bprecess.jl b/src/bprecess.jl
index d0f3be4..22b0592 100644
--- a/src/bprecess.jl
+++ b/src/bprecess.jl
@@ -177,8 +177,8 @@ http://adsabs.harvard.edu/abs/1983A%26A...128..263A.
 
 ### Example ###
 
-The SAO2000 catalogue gives the J2000 position and proper motion for the star HD
-119288.  Find the B1950 position.
+The SAO2000 catalogue gives the J2000 position and proper motion for the star
+HD 119288.  Find the B1950 position.
 
 * RA(2000) = 13h 42m 12.740s
 * Dec(2000) = 8d 23' 17.69''
@@ -210,4 +210,4 @@ and equinox of J2000.0" -- from the Explanatory Supplement (1992), p. 180
 
 Code of this function is based on IDL Astronomy User's Library.
 """
-bprecess
+function bprecess end
diff --git a/src/co_aberration.jl b/src/co_aberration.jl
index 4a498a3..2abc70a 100644
--- a/src/co_aberration.jl
+++ b/src/co_aberration.jl
@@ -63,8 +63,9 @@ julia> co_aberration(jd,ten(2,46,11.331)*15,ten(49,20,54.54))
 (30.04404628365077, 6.699400463119431)
 ```
 
-d_ra = 30.04404628365103'' (≈ 2.003s)
-d_dec = 6.699400463118504''
+d\\_ra = 30.04404628365103'' (≈ 2.003s)
+
+d\\_dec = 6.699400463118504''
 
 ### Notes ###
 
diff --git a/src/co_refract.jl b/src/co_refract.jl
index ccd9367..0c331a0 100644
--- a/src/co_refract.jl
+++ b/src/co_refract.jl
@@ -150,7 +150,7 @@ Code of this function is based on IDL Astronomy User's Library.
 function co_refract_forward(alt::Real, pre::Real, temp::Real)
 
     if alt < 15
-        ref =3.569*@evalpoly(alt, 0.1594, 0.0196, 0.00002)/@evalpoly(alt, 1, 0.505, 0.0845)
+        ref = 3.569*@evalpoly(alt, 0.1594, 0.0196, 0.00002)/@evalpoly(alt, 1, 0.505, 0.0845)
     else
         ref = 0.0166667 / tand((alt + 7.31 / (alt + 4.4)))
     end
diff --git a/src/ct2lst.jl b/src/ct2lst.jl
index 6485734..8425e7e 100644
--- a/src/ct2lst.jl
+++ b/src/ct2lst.jl
@@ -22,9 +22,9 @@ Convert from Local Civil Time to Local Mean Sidereal Time.
 The function can be called in two different ways.  The only argument common to
 both methods is `longitude`:
 
-* `longitude`: the longitude in degrees (east of Greenwich) of the place for which the local
-  sidereal time is desired.  The Greenwich mean sidereal time (GMST) can be found by setting
-  longitude = `0`.
+* `longitude`: the longitude in degrees (east of Greenwich) of the place for which
+  the local sidereal time is desired.  The Greenwich mean sidereal time (GMST) can
+  be found by setting longitude = `0`.
 
 The civil date to be converted to mean sidereal time can be specified either by
 providing the Julian days:
@@ -62,7 +62,7 @@ julia> lst = ct2lst(-76.72, -4, DateTime(2008, 7, 30, 15, 53))
 11.356505172312609
 
 julia> sixty(lst)
-3-element StaticArrays.SArray{Tuple{3},Float64,1,3} with indices SOneTo(3):
+3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
  11.0
  21.0
  23.418620325392112
@@ -76,17 +76,18 @@ Julian days.
 ```jldoctest
 julia> using AstroLib, Dates
 
-julia> longitude=ten(8, 43); # Convert longitude to decimals.
+julia> longitude = ten(8, 43); # Convert longitude to decimals.
+
+julia> # Get number of Julian days. Remember to subtract the time zone in
+       # order to convert local time to UTC.
 
 julia> jd = jdcnv(DateTime(2015, 11, 24, 13, 21) - Dates.Hour(1));
-# Get number of Julian days. Remember to subtract the time zone in
-# order to convert local time to UTC.
 
 julia> lst = ct2lst(longitude, jd) # Calculate Greenwich Mean Sidereal Time.
 17.140685171005316
 
 julia> sixty(lst)
-3-element StaticArrays.SArray{Tuple{3},Float64,1,3} with indices SOneTo(3):
+3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
  17.0
   8.0
  26.466615619137883
@@ -96,6 +97,8 @@ julia> sixty(lst)
 
 Code of this function is based on IDL Astronomy User's Library.
 """
+function ct2lst end
+
 ct2lst(long::Real, jd::Real) = ct2lst(promote(float(long), float(jd))...)
 
 function ct2lst(long::T, tz::T, date::DateTime) where {T<:AbstractFloat}
diff --git a/src/daycnv.jl b/src/daycnv.jl
index 87efaae..8a8f9a5 100644
--- a/src/daycnv.jl
+++ b/src/daycnv.jl
@@ -35,4 +35,4 @@ julia> daycnv(2440000)
 
 `jdcnv` is the inverse of this function.
 """
-const daycnv=Dates.julian2datetime
+const daycnv = Dates.julian2datetime
diff --git a/src/eci2geo.jl b/src/eci2geo.jl
index 44b7fd4..7f43f56 100644
--- a/src/eci2geo.jl
+++ b/src/eci2geo.jl
@@ -45,9 +45,9 @@ The three coordinates can be passed as a 3-tuple `(x, y, z)`.  In addition, `x`,
 
 The 3-tuple of geographical coordinate (latitude, longitude, altitude).
 
-* latitude: latitude, in degrees.
-* longitude: longitude, in degrees.
-* altitude: altitude, in kilometers.
+* `latitude`: latitude, in degrees.
+* `longitude`: longitude, in degrees.
+* `altitude`: altitude, in kilometers.
 
 If ECI coordinates are given as arrays, a 3-tuple of arrays of the same length
 is returned.
@@ -91,8 +91,7 @@ function eci2geo(x::AbstractArray{X}, y::AbstractArray{<:Real}, z::AbstractArray
     long = similar(x, typex)
     alt  = similar(x, typex)
     for i in eachindex(x)
-        lat[i], long[i], alt[i] =
-            eci2geo(x[i], y[i], z[i], jd[i])
+        lat[i], long[i], alt[i] = eci2geo(x[i], y[i], z[i], jd[i])
     end
     return lat, long, alt
 end
diff --git a/src/euler.jl b/src/euler.jl
index 2993ccf..fe5275d 100644
--- a/src/euler.jl
+++ b/src/euler.jl
@@ -104,7 +104,7 @@ function euler(ai::AbstractVector{R}, bi::AbstractVector{<:Real}, select::Intege
     bi_out = similar(bi, typeai)
     for i in eachindex(ai)
         ai_out[i], bi_out[i] = euler(ai[i], bi[i], select,
-                                        FK4=FK4, radians=radians)
+                                     FK4=FK4, radians=radians)
     end
     return ai_out, bi_out
 end
diff --git a/src/flux2mag.jl b/src/flux2mag.jl
index 8cc1dbd..af630fe 100644
--- a/src/flux2mag.jl
+++ b/src/flux2mag.jl
@@ -25,23 +25,25 @@ This is the reverse of `mag2flux`.
 * `flux`: the flux to be converted in magnitude, expressed in
   erg/(s cm² Å).
 * `zero_point`: the zero point level of the magnitude.  If not
- supplied then defaults to 21.1 (Code et al 1976).  Ignored if the `ABwave`
- keyword is supplied
+  supplied then defaults to 21.1 (Code et al 1976).  Ignored if the `ABwave`
+  keyword is supplied
 * `ABwave` (optional numeric keyword): wavelength in Angstroms.
- If supplied, then returns Oke AB magnitudes (Oke & Gunn 1983, ApJ, 266, 713;
- http://adsabs.harvard.edu/abs/1983ApJ...266..713O).
+  If supplied, then returns Oke AB magnitudes (Oke & Gunn 1983, ApJ, 266, 713;
+  http://adsabs.harvard.edu/abs/1983ApJ...266..713O).
 
 ### Output ###
 
 The magnitude.
 
 If the `ABwave` keyword is set then magnitude is given by the expression
-
-\$\$\\text{ABmag} = -2.5\\log_{10}(f) - 5\\log_{10}(\\text{ABwave}) - 2.406\$\$
+```math
+\\text{ABmag} = -2.5\\log_{10}(f) - 5\\log_{10}(\\text{ABwave}) - 2.406
+```
 
 Otherwise, magnitude is given by the expression
-
-\$\$\\text{mag} = -2.5\\log_{10}(\\text{flux}) - \\text{zero point}\$\$
+```math
+\\text{mag} = -2.5\\log_{10}(\\text{flux}) - \\text{zero point}
+```
 
 ### Example ###
 
diff --git a/src/gal_uvw.jl b/src/gal_uvw.jl
index 3ba5234..71fdab1 100644
--- a/src/gal_uvw.jl
+++ b/src/gal_uvw.jl
@@ -47,26 +47,26 @@ coordinates, (2) proper motion, (3) parallax, and (4) radial velocity.
 User must supply a position, proper motion, radial velocity and parallax.
 Either scalars or arrays all of the same length can be supplied.
 
-(1) Position:
+1. Position:
 
-* `ra`: right ascension, in degrees
-* `dec`: declination, in degrees
+   * `ra`: right ascension, in degrees
+   * `dec`: declination, in degrees
 
-(2) Proper Motion
+2. Proper Motion
 
-* `pmra`: proper motion in right ascension in arc units (typically
-  milli-arcseconds/yr).  If given \$\\mu_\\alpha\$ -- proper motion in seconds of
-  time/year -- then this is equal to \$15 \\mu_\\alpha \\cos(\\text{dec})\$.
-* `pmdec`: proper motion in declination (typically mas/yr).
+   * `pmra`: proper motion in right ascension in arc units (typically
+     milli-arcseconds/yr).  If given \$\\mu_\\alpha\$ -- proper motion in seconds of
+     time/year -- then this is equal to \$15 \\mu_\\alpha \\cos(\\text{dec})\$.
+   * `pmdec`: proper motion in declination (typically mas/yr).
 
-(3) Radial Velocity
+3. Radial Velocity
 
-* `vrad`: velocity in km/s
+   * `vrad`: velocity in km/s
 
-(4) Parallax
+4. Parallax
 
-* `plx`: parallax with same distance units as proper motion measurements
-  typically milliarcseconds (mas)
+   * `plx`: parallax with same distance units as proper motion measurements
+     typically milliarcseconds (mas)
 
 If you know the distance in parsecs, then set `plx` to \$1000/\\text{distance}\$,
 if proper motion measurements are given in milli-arcseconds/yr.
@@ -105,7 +105,7 @@ Hipparcos catalog, and correct to the LSR.
 ```jldoctest
 julia> using AstroLib
 
-julia> ra=ten(1,9,42.3)*15.; dec = ten(61,32,49.5);
+julia> ra = ten(1,9,42.3)*15.; dec = ten(61,32,49.5);
 
 julia> pmra = 627.89;  pmdec = 77.84; # mas/yr
 
diff --git a/src/gcirc.jl b/src/gcirc.jl
index d8fe3a3..f2922a5 100644
--- a/src/gcirc.jl
+++ b/src/gcirc.jl
@@ -86,7 +86,7 @@ julia> gcirc(0, 120, -43, 175, +22)
 ### Notes ###
 
 * The function `sphdist` provides an alternate method of computing a spherical
- distance.
+  distance.
 * The Haversine formula can give rounding errors for antipodal points.
 
 Code of this function is based on IDL Astronomy User's Library.
diff --git a/src/geo2eci.jl b/src/geo2eci.jl
index dc35ea4..a8da138 100644
--- a/src/geo2eci.jl
+++ b/src/geo2eci.jl
@@ -84,8 +84,7 @@ function geo2eci(lat::AbstractArray{LA}, long::AbstractArray{<:Real},
     y = similar(lat, typela)
     z = similar(lat, typela)
     for i in eachindex(lat)
-        x[i], y[i], z[i] =
-            geo2eci(lat[i], long[i], alt[i], jd[i])
+        x[i], y[i], z[i] = geo2eci(lat[i], long[i], alt[i], jd[i])
     end
     return x, y, z
 end
diff --git a/src/geo2geodetic.jl b/src/geo2geodetic.jl
index 13a6549..16bd17c 100644
--- a/src/geo2geodetic.jl
+++ b/src/geo2geodetic.jl
@@ -141,7 +141,7 @@ planetodetic) to geographic coordinates, can be used to estimate the accuracy of
 julia> using AstroLib
 
 julia> collect(geodetic2geo(geo2geodetic(67.2, 13.4, 1.2))) - [67.2, 13.4, 1.2]
-3-element Array{Float64,1}:
+3-element Vector{Float64}:
  -3.5672513831741526e-9
   0.0
   9.484211194177306e-10
diff --git a/src/geo2mag.jl b/src/geo2mag.jl
index ea05a02..4e55112 100644
--- a/src/geo2mag.jl
+++ b/src/geo2mag.jl
@@ -70,7 +70,7 @@ geomagnetic coordinates in 2016:
 julia> using AstroLib
 
 julia> geo2mag(ten(35,0,42), ten(135,46,6), 2016)
-(36.86579228937769, -60.184060536651614)
+(36.86579228937769, -60.1840605366516)
 ```
 
 ### Notes ###
diff --git a/src/get_date.jl b/src/get_date.jl
index ad1d57e..c8955c7 100644
--- a/src/get_date.jl
+++ b/src/get_date.jl
@@ -55,11 +55,11 @@ julia> get_date(DateTime(21937, 05, 30, 09, 59, 00), timetag=true)
 ### Notes ###
 
 1. A discussion of the DATExxx syntax in FITS headers can be found in
- http://www.cv.nrao.edu/fits/documents/standards/year2000.txt
+   http://www.cv.nrao.edu/fits/documents/standards/year2000.txt
 
 2. Those who wish to use need further flexibility in their date formats (e.g. to
- use TAI time) should look at Bill Thompson's time routines in
- http://sohowww.nascom.nasa.gov/solarsoft/gen/idl/time
+   use TAI time) should look at Bill Thompson's time routines in
+   http://sohowww.nascom.nasa.gov/solarsoft/gen/idl/time
 """
 function get_date(dt::DateTime, old::Bool, timetag::Bool)
     # Based on `Base.string' definition in base/dates/io.jl.
diff --git a/src/helio.jl b/src/helio.jl
index 023cb03..ba62a36 100644
--- a/src/helio.jl
+++ b/src/helio.jl
@@ -13,7 +13,7 @@ const dpd = @SMatrix [ 0.00000037  0.00001906 -0.00594749 -0.12534081  0.1604768
                       -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")
+                    6=>"saturn", 7=>"uranus", 8=>"neptune", 9=>"pluto")
 
 function _helio(jd::T, num::Integer, radians::Bool) where {T<:AbstractFloat}
 
@@ -63,8 +63,8 @@ for Saturn.
 * `jd`: julian date, scalar or vector
 * `num`: integer denoting planet number, scalar or vector
   1 = Mercury, 2 = Venus, ... 9 = Pluto
-* `radians`(optional): if this keyword is set to
-  `true`, than the longitude and latitude output are in radians rather than degrees.
+* `radians`(optional): if this keyword is set to `true`, then
+  the longitude and latitude output are in radians rather than degrees.
 
 ### Output ###
 
@@ -74,32 +74,32 @@ for Saturn.
 
 ### Example ###
 
-(1) Find heliocentric position of Venus on August 23, 2000
+- Find heliocentric position of Venus on August 23, 2000
 
-```jldoctest
-julia> using AstroLib
+  ```jldoctest
+  julia> using AstroLib
 
-julia> helio(jdcnv(2000,08,23,0), 2)
-(0.7213758288364316, 198.39093251916148, 2.887355631705488)
-```
+  julia> helio(jdcnv(2000,08,23,0), 2)
+  (0.7213758288364316, 198.39093251916148, 2.887355631705488)
+  ```
 
-(2) Find the current heliocentric positions of all the planets
+- Find the current heliocentric positions of all the planets
 
-```jldoctest
-julia> using AstroLib
+  ```jldoctest
+  julia> using AstroLib
 
-julia> helio.([jdcnv(1900)], 1:9)
-9-element Array{Tuple{Float64,Float64,Float64},1}:
- (0.4207394142180803, 202.60972662618906, 3.0503005607270532)
- (0.7274605731764012, 344.5381482401048, -3.3924346961624785)
- (0.9832446886519147, 101.54969268801035, 0.012669354526696368)
- (1.4212659241051142, 287.8531100442217, -1.5754626002228043)
- (5.386813769590955, 235.91306092135062, 0.9131692817310215)
- (10.054339927304339, 268.04069870870387, 1.0851704598594278)
- (18.984683376211326, 250.0555468087738, 0.05297087029604253)
- (29.87722677219009, 87.07244903504716, -1.245060583142733)
- (46.9647515992327, 75.94692594417324, -9.576681044165511)
-```
+  julia> helio.([jdcnv(1900)], 1:9)
+  9-element Vector{Tuple{Float64, Float64, Float64}}:
+   (0.4207394142180803, 202.60972662618906, 3.0503005607270532)
+   (0.7274605731764012, 344.5381482401048, -3.3924346961624785)
+   (0.9832446886519147, 101.54969268801035, 0.012669354526696368)
+   (1.4212659241051142, 287.8531100442217, -1.5754626002228043)
+   (5.386813769590955, 235.91306092135062, 0.9131692817310215)
+   (10.054339927304339, 268.04069870870387, 1.0851704598594278)
+   (18.984683376211326, 250.0555468087738, 0.05297087029604253)
+   (29.87722677219009, 87.07244903504716, -1.245060583142733)
+   (46.9647515992327, 75.94692594417324, -9.576681044165511)
+  ```
 ### Notes ###
 
 This program is based on the two-body model and thus neglects
@@ -122,8 +122,7 @@ function helio(jd::AbstractVector{P}, num::AbstractVector{<:Real},
     hlong_out = similar(jd,  typejd)
     hlat_out = similar(jd,  typejd)
     for i in eachindex(jd)
-        hrad_out[i], hlong_out[i], hlat_out[i] =
-        helio(jd[i], num[i], radians)
+        hrad_out[i], hlong_out[i], hlat_out[i] = helio(jd[i], num[i], radians)
     end
     return hrad_out, hlong_out, hlat_out
 end
diff --git a/src/helio_rv.jl b/src/helio_rv.jl
index 071eba1..83a0ee1 100644
--- a/src/helio_rv.jl
+++ b/src/helio_rv.jl
@@ -38,33 +38,33 @@ date(s) specified by `jd`.
 
 ### Example ###
 
-(1) What was the heliocentric radial velocity of the primary component of HU Tau
-at 1730 UT 25 Oct 1994?
+- What was the heliocentric radial velocity of the primary component of HU Tau
+  at 1730 UT 25 Oct 1994?
 
-```jldoctest
-julia> using AstroLib
+  ```jldoctest
+  julia> using AstroLib
 
-julia> jd = juldate(94, 10, 25, 17, 30); # Obtain Geocentric Julian days
+  julia> jd = juldate(94, 10, 25, 17, 30); # Obtain Geocentric Julian days
 
-julia> hjd = helio_jd(jd, ten(04, 38, 16) * 15, ten(20, 41, 05)); # Convert to HJD
+  julia> hjd = helio_jd(jd, ten(04, 38, 16) * 15, ten(20, 41, 05)); # Convert to HJD
 
-julia> helio_rv(hjd, 46487.5303, 2.0563056, -6, 59.3)
--62.965570107789475
-```
+  julia> helio_rv(hjd, 46487.5303, 2.0563056, -6, 59.3)
+  -62.965570107789475
+  ```
 
-NB: the functions `juldate` and `helio_jd` return a reduced HJD (HJD - 2400000)
-and so T and P must be specified in the same fashion.
+  NB: the functions `juldate` and `helio_jd` return a reduced HJD (HJD - 2400000)
+  and so T and P must be specified in the same fashion.
 
-(2) Plot two cycles of an eccentric orbit, \$e=0.6\$, \$\\omega=45\\degree\$ for
-both components of a binary star.  Use
-[PyPlot.jl](https://github.com/JuliaPy/PyPlot.jl) for plotting.
+- Plot two cycles of an eccentric orbit, \$e=0.6\$, \$\\omega=45\\degree\$ for
+  both components of a binary star.  Use
+  [PyPlot.jl](https://github.com/JuliaPy/PyPlot.jl) for plotting.
 
-```julia
-using PyPlot
-φ = range(0, stop=2, length=1000); # Generate 1000 phase points
-plot(φ ,helio_rv.(φ, 0, 1, 0, 100, 0.6, 45)) # Plot 1st component
-plot(φ ,helio_rv.(φ, 0, 1, 0, 100, 0.6, 45+180)) # Plot 2nd component
-```
+  ```julia
+  using PyPlot
+  φ = range(0, stop=2, length=1000); # Generate 1000 phase points
+  plot(φ, helio_rv.(φ, 0, 1, 0, 100, 0.6, 45)) # Plot 1st component
+  plot(φ, helio_rv.(φ, 0, 1, 0, 100, 0.6, 45+180)) # Plot 2nd component
+  ```
 
 ### Notes ###
 
diff --git a/src/hor2eq.jl b/src/hor2eq.jl
index 02d8e0d..ec81e13 100644
--- a/src/hor2eq.jl
+++ b/src/hor2eq.jl
@@ -54,13 +54,16 @@ hor2eq(alt::Real, az::Real, jd::Real, obsname::AbstractString; kwargs...) =
 
 
 """
-    hor2eq(alt, az, jd[, obsname; ws=false, B1950=false, precession=true, nutate=true,
-           aberration=true, refract=true, lat=NaN, lon=NaN, altitude=0, pressure=NaN,
-           temperature=NaN]) -> ra, dec, ha
+    hor2eq(alt, az, jd[, obsname;
+           ws=false, B1950=false, precession=true, nutate=true,
+           aberration=true, refract=true,
+           lat=NaN, lon=NaN, altitude=0,
+           pressure=NaN, temperature=NaN]) -> ra, dec, ha
 
-    hor2eq(alt, az, jd, lat, lon[, altitude=0; ws=false, B1950=false,
-           precession=true, nutate=true, aberration=true, refract=true, pressure=NaN,
-           temperature=NaN]) -> ra, dec, ha
+    hor2eq(alt, az, jd, lat, lon[, altitude=0;
+           ws=false, B1950=false, precession=true, nutate=true,
+           aberration=true, refract=true,
+           pressure=NaN, temperature=NaN]) -> ra, dec, ha
 
 
 ### Purpose
@@ -135,7 +138,7 @@ julia> using AstroLib
 
 julia> ra_o, dec_o = hor2eq(ten(37,54,41), ten(264,55,06), 2466879.7083333,
                             "kpno", pressure = 781, temperature = 273)
-(3.3224480269254717, 15.19061543702944, 54.61174536229464)
+(3.3224480269254713, 15.19061543702944, 54.61174536229464)
 
 julia> adstring(ra_o, dec_o)
 " 00 13 17.4  +15 11 26"
@@ -145,7 +148,7 @@ julia> adstring(ra_o, dec_o)
 
 Code of this function is based on IDL Astronomy User's Library.
 """
-hor2eq
+function hor2eq end
 
 # TODO: Make hor2eq type-stable, which it isn't currently because of keyword arguments
 # Note that the inner function `_hor2eq` is type stable
diff --git a/src/imf.jl b/src/imf.jl
index 9c043f2..12b2872 100644
--- a/src/imf.jl
+++ b/src/imf.jl
@@ -66,7 +66,7 @@ Show the number of stars per unit mass interval at 3 Msun for a Salpeter
 julia> using AstroLib
 
 julia> imf([3], [-1.35], [0.1, 110]) / 3
-1-element Array{Float64,1}:
+1-element Vector{Float64}:
  0.01294143518151214
 ```
 
@@ -74,5 +74,7 @@ julia> imf([3], [-1.35], [0.1, 110]) / 3
 
 Code of this function is based on IDL Astronomy User's Library.
 """
-imf(mass::AbstractVector{<:Real}, expon::AbstractVector{<:Real}, mass_range::AbstractVector{<:Real}) =
-    imf(float(mass), float(expon), float(mass_range))
+imf(mass::AbstractVector{<:Real},
+    expon::AbstractVector{<:Real},
+    mass_range::AbstractVector{<:Real}) =
+        imf(float(mass), float(expon), float(mass_range))
diff --git a/src/ismeuv.jl b/src/ismeuv.jl
index 61b3f6c..c991ee9 100644
--- a/src/ismeuv.jl
+++ b/src/ismeuv.jl
@@ -8,7 +8,7 @@ function ismeuv(wave::T, hcol::T, he1col::T, he2col::T, fano::Bool) where {T<:Ab
     # minexp = -708.39642  from alog((machar(double=double).xmin)) where double = 1b in IDL
     z = sqrt(r/(1 - r))
     denom = -expm1(-2*T(π)*z)
-    tauH = ((3.44e-16)*(r^4)*exp(-4*z*atan(1/z))*hcol)/denom
+    tauH = (3.44e-16 * r^4 *exp(-4*z*atan(1/z)) * hcol) / denom
     r *= 4
     tauHe2 = zero(T)
     if r < 1
@@ -78,12 +78,11 @@ The EUV optical depth is computed from the photoionization of hydrogen and heliu
 One has a model EUV spectrum with wavelength, w (in Angstroms).
 Find the EUV optical depth by 1e18 cm-2 of HI, with N(HeI)/N(HI) = N(HeII)/N(HI) = 0.05.
 
-
 ```jldoctest
 julia> using AstroLib
 
 julia> ismeuv.([670, 910], 1e19, 5e17, 5e17)
-2-element Array{Float64,1}:
+2-element Vector{Float64}:
  27.35393320556168
  62.683796028917286
 ```
diff --git a/src/kepler_solver.jl b/src/kepler_solver.jl
index c641704..5e434df 100644
--- a/src/kepler_solver.jl
+++ b/src/kepler_solver.jl
@@ -51,7 +51,9 @@ In order to find the position of a body in elliptic motion (e.g., in the
 two-body problem) at a given time \$t\$, one has to solve the [Kepler's
 equation](https://en.wikipedia.org/wiki/Kepler%27s_equation)
 
-``M(t) = E(t) - e\\sin E(t)``
+```math
+M(t) = E(t) - e \\sin E(t)
+```
 
 where \$M(t) = (t - t_{0})/P\$ is the mean anomaly, \$E(t)\$ the eccentric
 anomaly, \$e\$ the eccentricity of the orbit, \$t_0\$ is the time of periapsis
@@ -66,7 +68,7 @@ given and one wants to find the eccentric anomaly \$E(t)\$ at a specific time
 
 ### Output ###
 
-The eccentric anomaly \$E\$, restricted to the range \$[-\\pi, \\pi]\$.
+The eccentric anomaly \$E\$, restricted to the range \$[-π, π]\$.
 
 ### Method ###
 
@@ -80,31 +82,33 @@ motion \$0 \\leq e \\leq 1\$.
 
 ### Example ###
 
-(1) Find the eccentric anomaly for an orbit with eccentricity \$e = 0.7\$ and
-for \$M(t) = 8\\pi/3\$.
+- Find the eccentric anomaly for an orbit with eccentricity \$e = 0.7\$ and
+  for \$M(t) = 8π/3\$.
 
-```jldoctest
-julia> using AstroLib
+  ```jldoctest
+  julia> using AstroLib
 
-julia> ecc = 0.7;
+  julia> ecc = 0.7;
 
-julia> E = kepler_solver(8pi/3, ecc)
-2.5085279492864223
-```
+  julia> E = kepler_solver(8pi/3, ecc)
+  2.5085279492864223
+  ```
 
-(2) Plot the eccentric anomaly as a function of mean anomaly for eccentricity
-\$e = 0\$, \$0.5\$, \$0.9\$.  Recall that `kepler_solver` gives \$E \\in [-\\pi,
-\\pi]\$, use `mod2pi` to have it in \$[0, 2\\pi]\$.  Use
-[PyPlot.jl](https://github.com/JuliaPlots/Plots.jl/) for plotting.
+- Plot the eccentric anomaly as a function of mean anomaly for eccentricity
+  \$e = 0\$, \$0.5\$, \$0.9\$.  Recall that `kepler_solver` gives \$E ∈ [-π,
+  π]\$, use `mod2pi` to have it in \$[0, 2\\pi]\$.  Use
+  [PyPlot.jl](https://github.com/JuliaPlots/Plots.jl/) for plotting.
 
-```julia
-using AstroLib, PyPlot
-M = range(0, stop=2pi, length=1001)[1:end-1];
-for ecc in (0, 0.5, 0.9); plot(M, mod2pi.(kepler_solver.(M, ecc))); end
-```
+  ```julia
+  using AstroLib, PyPlot
+  M = range(0, stop=2pi, length=1001)[1:end-1];
+  for ecc in (0, 0.5, 0.9)
+      plot(M, mod2pi.(kepler_solver.(M, ecc)))
+  end
+  ```
 
 ### Notes ###
 
 The true anomaly can be calculated with `trueanom` function.
 """
-kepler_solver
+function kepler_solver end
diff --git a/src/mag2flux.jl b/src/mag2flux.jl
index 462defb..8332843 100644
--- a/src/mag2flux.jl
+++ b/src/mag2flux.jl
@@ -24,10 +24,10 @@ This is the reverse of `flux2mag`.
 
 * `mag`: the magnitude to be converted in flux.
 * `zero_point`: the zero point level of the magnitude.  If not supplied then defaults to
- 21.1 (Code et al 1976).  Ignored if the `ABwave` keyword is supplied
+  21.1 (Code et al 1976).  Ignored if the `ABwave` keyword is supplied
 * `ABwave` (optional numeric keyword): wavelength, in Angstroms.  If supplied, then the
- input `mag` is assumed to contain Oke AB magnitudes (Oke & Gunn 1983, ApJ, 266, 713;
- http://adsabs.harvard.edu/abs/1983ApJ...266..713O).
+  input `mag` is assumed to contain Oke AB magnitudes (Oke & Gunn 1983, ApJ, 266, 713;
+  http://adsabs.harvard.edu/abs/1983ApJ...266..713O).
 
 ### Output ###
 
@@ -53,7 +53,7 @@ julia> mag2flux(8.3, 12)
 7.58577575029182e-9
 
 julia> mag2flux(8.3, ABwave=12)
-3.6244115683017193e-7
+3.624411568301719e-7
 ```
 
 ### Notes ###
diff --git a/src/month_cnv.jl b/src/month_cnv.jl
index 0e112f0..38cccd7 100644
--- a/src/month_cnv.jl
+++ b/src/month_cnv.jl
@@ -43,13 +43,13 @@ format of the month name is influenced by the optional keywords.
 julia> using AstroLib
 
 julia> month_cnv.(["janua", "SEP", "aUgUsT"])
-3-element Array{Int64,1}:
+3-element Vector{Int64}:
  1
  9
  8
 
 julia> month_cnv.([2, 12, 6], short=true, low=true)
-3-element Array{String,1}:
+3-element Vector{String}:
  "feb"
  "dec"
  "jun"
diff --git a/src/moonpos.jl b/src/moonpos.jl
index 2ab16e0..ed51184 100644
--- a/src/moonpos.jl
+++ b/src/moonpos.jl
@@ -175,32 +175,34 @@ difference in the position calculation.
 
 ### Example ###
 
-(1) Find the position of the moon on April 12, 1992
+- Find the position of the moon on April 12, 1992
 
-```jldoctest
-julia> using AstroLib
+  ```jldoctest
+  julia> using AstroLib
 
-julia> jd = jdcnv(1992, 4, 12);
+  julia> jd = jdcnv(1992, 4, 12);
 
-julia> adstring(moonpos(jd)[1:2],precision=1)
-" 08 58 45.23  +13 46 06.1"
-```
+  julia> adstring(moonpos(jd)[1:2],precision=1)
+  " 08 58 45.23  +13 46 06.1"
+  ```
 
-This is within 1" from the position given in the Astronomical Almanac.
+  This is within 1" from the position given in the Astronomical Almanac.
 
-(2) Plot the Earth-moon distance during 2016 with sampling of 6 hours.  Use
-[PyPlot.jl](https://github.com/JuliaPlots/Plots.jl/) for plotting
+- Plot the Earth-moon distance during 2016 with sampling of 6 hours.  Use
+  [PyPlot.jl](https://github.com/JuliaPlots/Plots.jl/) for plotting
 
-```julia
-using PyPlot
-points = DateTime(2016):Dates.Hour(6):DateTime(2017);
-plot(points, moonpos(jdcnv.(points))[3])
-```
+  ```julia
+  using PyPlot
+  points = DateTime(2016):Dates.Hour(6):DateTime(2017);
+  plot(points, moonpos(jdcnv.(points))[3])
+  ```
 
 ### Notes ###
 
 Code of this function is based on IDL Astronomy User's Library.
 """
+function moonpos end
+
 moonpos(jd::Real; radians::Bool=false) = _moonpos(float(jd), radians)
 
 function moonpos(jd::AbstractArray{J}; radians::Bool=false) where {J<:Real}
diff --git a/src/nutate.jl b/src/nutate.jl
index 4599de0..027240a 100644
--- a/src/nutate.jl
+++ b/src/nutate.jl
@@ -84,33 +84,33 @@ If `jd` is an array, `long` and `obl` are arrays of the same length.
 
 ### Method ###
 
-Uses the formula in Chapter 22 of ``Astronomical Algorithms'' by Jean Meeus
+Uses the formula in Chapter 22 of "Astronomical Algorithms" by Jean Meeus
 (1998, 2nd ed.) which is based on the 1980 IAU Theory of Nutation and includes
 all terms larger than 0.0003".
 
 ### Example ###
 
-(1) Find the nutation in longitude and obliquity 1987 on Apr 10 at Oh.  This is
-example 22.a from Meeus
+- Find the nutation in longitude and obliquity 1987 on Apr 10 at 0h.  This is
+  example 22.a from Meeus
 
-```jldoctest
-julia> using AstroLib
+  ```jldoctest
+  julia> using AstroLib
 
-julia> jd = jdcnv(1987, 4, 10);
+  julia> jd = jdcnv(1987, 4, 10);
 
-julia> nutate(jd)
-(-3.787931077110494, 9.44252069864449)
-```
+  julia> nutate(jd)
+  (-3.787931077110494, 9.44252069864449)
+  ```
 
-(2) Plot the daily nutation in longitude and obliquity during the 21st century.
-Use [PyPlot.jl](https://github.com/JuliaPlots/Plots.jl/) for plotting.
+- Plot the daily nutation in longitude and obliquity during the 21st century.
+  Use [PyPlot.jl](https://github.com/JuliaPlots/Plots.jl/) for plotting.
 
-```julia
-using PyPlot
-years = DateTime(2000):DateTime(2100);
-long, obl = nutate(jdcnv.(years));
-plot(years, long); plot(years, obl)
-```
+  ```julia
+  using PyPlot
+  years = DateTime(2000):DateTime(2100);
+  long, obl = nutate(jdcnv.(years));
+  plot(years, long); plot(years, obl)
+  ```
 
 You can see both the dominant large scale period of nutation, of 18.6 years, and
 smaller oscillations with shorter periods.
diff --git a/src/ordinal.jl b/src/ordinal.jl
index eec4f6a..8ec764a 100644
--- a/src/ordinal.jl
+++ b/src/ordinal.jl
@@ -8,15 +8,15 @@
 Convert an integer to a correct English ordinal string.
 
 ### Explanation ###
-The first four ordinal strings are "1st", "2nd", "3rd", "4th" ....
+The first four ordinal strings are "1st", "2nd", "3rd", "4th", ....
 
 ### Arguments ###
 
-* `num`: number to be made ordinal. It should be of type int.
+* `num`: number to be made ordinal. It should be of type `Integer`.
 
 ### Output ###
 
-* `result`: ordinal string, such as '1st' '3rd '164th' '87th' etc
+* `result`: ordinal string, such as "1st", "3rd", "164th", "87th", etc.
 
 ### Example ###
 
@@ -24,7 +24,7 @@ The first four ordinal strings are "1st", "2nd", "3rd", "4th" ....
 julia> using AstroLib
 
 julia> ordinal.(1:5)
-5-element Array{String,1}:
+5-element Vector{String}:
  "1st"
  "2nd"
  "3rd"
diff --git a/src/paczynski.jl b/src/paczynski.jl
index 027af48..70c5230 100644
--- a/src/paczynski.jl
+++ b/src/paczynski.jl
@@ -70,7 +70,7 @@ Calculate the microlensing amplification for \$u = 10^{-10}\$, \$10^{-1}\$,
 julia> using AstroLib
 
 julia> paczynski.([1e-10, 1e-1, 1, 10, 1e10])
-5-element Array{Float64,1}:
+5-element Vector{Float64}:
   1.0e10
  10.037461005722337
   1.3416407864998738
diff --git a/src/planck_freq.jl b/src/planck_freq.jl
index f12e54e..4c0ee6d 100644
--- a/src/planck_freq.jl
+++ b/src/planck_freq.jl
@@ -19,7 +19,9 @@ Calculate the flux of a black body per unit frequency.
 Return the spectral radiance of a black body per unit frequency using [Planck's
 law](https://en.wikipedia.org/wiki/Planck%27s_law)
 
-`` B_\\nu(\\nu, T) = \\frac{2h\\nu ^3}{c^2} \\frac{1}{e^\\frac{h\\nu}{k_\\mathrm{B}T} - 1} ``
+```math
+B_\\nu(\\nu, T) = \\frac{2h\\nu ^3}{c^2} \\frac{1}{e^\\frac{h\\nu}{k_\\mathrm{B}T} - 1}
+```
 
 ### Arguments ###
 
diff --git a/src/planck_wave.jl b/src/planck_wave.jl
index 56cf9a5..9363ad5 100644
--- a/src/planck_wave.jl
+++ b/src/planck_wave.jl
@@ -19,7 +19,9 @@ Calculate the flux of a black body per unit wavelength.
 Return the spectral radiance of a black body per unit wavelength using [Planck's
 law](https://en.wikipedia.org/wiki/Planck%27s_law)
 
-`` B_\\lambda(\\lambda, T) =\\frac{2hc^2}{\\lambda^5}\\frac{1}{e^{\\frac{hc}{\\lambda k_\\mathrm{B}T}} - 1} ``
+```math
+B_λ(λ, T) = \\frac{2hc^2}{λ^5}\\frac{1}{e^{\\frac{hc}{λ k_\\mathrm{B}T}} - 1}
+```
 
 ### Arguments ###
 
diff --git a/src/posang.jl b/src/posang.jl
index 77c47d9..1bef56f 100644
--- a/src/posang.jl
+++ b/src/posang.jl
@@ -51,14 +51,14 @@ center.
 
 ### Arguments ###
 
-* `units`: integer, can be either 0, or 1, or 2.  Describes units of inputs and
- output:
-    * 0: everything (input right ascensions and declinations, and output
-      distance) is radians
-    * 1: right ascensions are in decimal hours, declinations in decimal degrees,
-      output distance in degrees
-    * 2: right ascensions and declinations are in degrees, output distance in
-      degrees
+* `units`: integer, can be either 0, or 1, or 2.  Describes units of inputs
+  and output:
+  * 0: everything (input right ascensions and declinations, and output
+    distance) is radians
+  * 1: right ascensions are in decimal hours, declinations in decimal degrees,
+    output distance in degrees
+  * 2: right ascensions and declinations are in degrees, output distance in
+    degrees
 * `ra1`:  right ascension or longitude of point 1
 * `dec1`: declination or latitude of point 1
 * `ra2`: right ascension or longitude of point 2
diff --git a/src/precess_cd.jl b/src/precess_cd.jl
index 7846769..c89e4d3 100644
--- a/src/precess_cd.jl
+++ b/src/precess_cd.jl
@@ -41,7 +41,7 @@ The coordinate matrix is precessed from epoch1 to epoch2.
 
 ### Arguments ###
 
-* `cd`: 2 x 2 coordinate description matrix in degrees
+* `cd`: 2×2 coordinate description matrix in degrees
 * `epoch1`: original equinox of coordinates, scalar
 * `epoch2`: equinox of precessed coordinates, scalar
 * `crval_old`: 2 element vector containing right ascension and declination
@@ -61,7 +61,7 @@ The coordinate matrix is precessed from epoch1 to epoch2.
 julia> using AstroLib
 
 julia> precess_cd([20 60; 45 45], 1950, 2000, [34, 58], [12, 83])
-2×2 Array{Float64,2}:
+2×2 Matrix{Float64}:
   48.8944  147.075
  110.188   110.365
 ```
diff --git a/src/premat.jl b/src/premat.jl
index 002a3e7..4e04fd2 100644
--- a/src/premat.jl
+++ b/src/premat.jl
@@ -59,8 +59,8 @@ Return the precession matrix from 1950.0 to 1975.0 in the FK4 system
 ```jldoctest
 julia> using AstroLib
 
-julia> premat(1950,1975,FK4=true)
-3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9} with indices SOneTo(3)×SOneTo(3):
+julia> premat(1950, 1975, FK4=true)
+3×3 StaticArraysCore.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
  0.999981    -0.00558775  -0.00242909
  0.00558775   0.999984    -6.78691e-6
  0.00242909  -6.78633e-6   0.999997
diff --git a/src/sixty.jl b/src/sixty.jl
index 774804e..5e6a48f 100644
--- a/src/sixty.jl
+++ b/src/sixty.jl
@@ -38,7 +38,7 @@ An array of three `AbstractFloat`, that are the sexagesimal counterpart
 julia> using AstroLib
 
 julia> sixty(-0.615)
-3-element StaticArrays.SArray{Tuple{3},Float64,1,3} with indices SOneTo(3):
+3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
  -0.0
  36.0
  54.0
diff --git a/src/sunpos.jl b/src/sunpos.jl
index fd6ad54..d4c691b 100644
--- a/src/sunpos.jl
+++ b/src/sunpos.jl
@@ -62,8 +62,20 @@ function sunpos(jd::Real; radians::Bool = false)
     end
 end
 
+function sunpos(jd::AbstractVector{J}; radians::Bool=false) where {J<:Real}
+    typej = float(J)
+    ra = similar(jd, typej)
+    dec = similar(jd, typej)
+    longmed = similar(jd, typej)
+    oblt = similar(jd, typej)
+    for i in eachindex(jd)
+        ra[i], dec[i], longmed[i], oblt[i] = sunpos(jd[i], radians=radians)
+    end
+    return ra, dec, longmed, oblt
+end
+
 """
-    sunpos(jd[, radians=true]) -> ra, dec, elong, obliquity
+    sunpos(jd[, radians=false]) -> ra, dec, elong, obliquity
 
 ### Purpose ###
 
@@ -98,30 +110,30 @@ by CD Pike, which was adapted from a FORTRAN routine by B. Emerson (RGO).
 
 ### Example ###
 
-(1) Find the apparent right ascension and declination of the Sun on May 1, 1982
+- Find the apparent right ascension and declination of the Sun on May 1, 1982
 
-```jldoctest
-julia> using AstroLib
+  ```jldoctest
+  julia> using AstroLib
 
-julia> adstring(sunpos(jdcnv(1982, 5, 1))[1:2], precision=2)
-" 02 31 32.614  +14 54 34.92"
-```
+  julia> adstring(sunpos(jdcnv(1982, 5, 1))[1:2], precision=2)
+  " 02 31 32.614  +14 54 34.92"
+  ```
 
-The Astronomical Almanac gives `02 31 32.58 +14 54 34.9` so the error for this
-case is < 0.5".
+  The Astronomical Almanac gives `02 31 32.58 +14 54 34.9` so the error for this
+  case is < 0.5".
 
-(2) Plot the apparent right ascension, in hours, and declination of the Sun, in
-degrees, for every day in 2016.  Use
-[PyPlot.jl](https://github.com/JuliaPlots/Plots.jl/) for plotting.
+- Plot the apparent right ascension, in hours, and declination of the Sun, in
+  degrees, for every day in 2016.  Use
+  [PyPlot.jl](https://github.com/JuliaPlots/Plots.jl/) for plotting.
 
-```julia
-using PyPlot
-using Dates
+  ```julia
+  using PyPlot
+  using Dates
 
-days = DateTime(2016):Day(1):DateTime(2016, 12, 31);
-ra, declin = sunpos(jdcnv.(days));
-plot(days, ra/15); plot(days, declin)
-```
+  days = DateTime(2016):Day(1):DateTime(2016, 12, 31);
+  ra, declin = sunpos(jdcnv.(days));
+  plot(days, ra/15); plot(days, declin)
+  ```
 
 ### Notes ###
 
@@ -134,14 +146,4 @@ The returned `ra` and `dec` are in the given date's equinox.
 
 Code of this function is based on IDL Astronomy User's Library.
 """
-function sunpos(jd::AbstractVector{J}; radians::Bool=false) where {J<:Real}
-    typej = float(J)
-    ra = similar(jd, typej)
-    dec = similar(jd, typej)
-    longmed = similar(jd, typej)
-    oblt = similar(jd, typej)
-    for i in eachindex(jd)
-        ra[i], dec[i], longmed[i], oblt[i] = sunpos(jd[i], radians=radians)
-    end
-    return ra, dec, longmed, oblt
-end
+function sunpos end
diff --git a/src/ten.jl b/src/ten.jl
index 2c303a4..cba874c 100644
--- a/src/ten.jl
+++ b/src/ten.jl
@@ -80,4 +80,4 @@ If it is important to give sense to negative zero, you can either make sure to
 pass a floating point negative zero `-0.0` (this is the best option), or use
 negative minutes and seconds, or non-integer negative degrees and minutes.
 """
-ten
+function ten end
diff --git a/src/trueanom.jl b/src/trueanom.jl
index dd166f5..c614e8d 100644
--- a/src/trueanom.jl
+++ b/src/trueanom.jl
@@ -18,15 +18,19 @@ Calculate true anomaly for a particle in elliptic orbit with eccentric anomaly
 
 In the two-body problem, once that the [Kepler's
 equation](https://en.wikipedia.org/wiki/Kepler%27s_equation) is solved and
-\$E(t)\$ is determined, the polar coordinates \$(r(t), \\theta(t))\$ of the body
+\$E(t)\$ is determined, the polar coordinates \$(r(t), θ(t))\$ of the body
 at time \$t\$ in the elliptic orbit are given by
 
-`` \\theta(t) = 2\\arctan \\left(\\sqrt{\\frac{1 + e}{1 - e}} \\tan\\frac{E(t)}{2} \\right)``
+```math
+θ(t) = 2\\arctan \\left(\\sqrt{\\frac{1 + e}{1 - e}} \\tan\\frac{E(t)}{2} \\right)
+```
 
-`` r(t) = \\frac{a(1 - e^{2})}{1 + e\\cos(\\theta(t) - \\theta_{0})}``
+```math
+r(t) = \\frac{a(1 - e^{2})}{1 + e\\cos(θ(t) - θ_0)}
+```
 
-in which \$a\$ is the semi-major axis of the orbit, and \$\\theta_0\$ the value
-of angular coordinate at time \$t = t_{0}\$.
+in which \$a\$ is the semi-major axis of the orbit, and \$θ_0\$ the value
+of angular coordinate at time \$t = t_0\$.
 
 ### Arguments ###
 
diff --git a/src/uvbybeta.jl b/src/uvbybeta.jl
index dace68a..6ed9d0e 100644
--- a/src/uvbybeta.jl
+++ b/src/uvbybeta.jl
@@ -255,14 +255,14 @@ Derive dereddened colors, metallicity, and Teff from Stromgren colors.
 * `c1`: Stromgren Balmer discontinuity parameter
 * `n`: Integer which can be any value between 1 to 8, giving approximate stellar
   classification.
-  (1) B0 - A0, classes III - V, 2.59 < Hbeta < 2.88,-0.20 <   c0   < 1.00
-  (2) B0 - A0, class   Ia     , 2.52 < Hbeta < 2.59,-0.15 <   c0   < 0.40
-  (3) B0 - A0, class   Ib     , 2.56 < Hbeta < 2.61,-0.10 <   c0   < 0.50
-  (4) B0 - A0, class   II     , 2.58 < Hbeta < 2.63,-0.10 <   c0   < 0.10
-  (5) A0 - A3, classes III - V, 2.87 < Hbeta < 2.93,-0.01 < (b-y)o < 0.06
-  (6) A3 - F0, classes III - V, 2.72 < Hbeta < 2.88, 0.05 < (b-y)o < 0.22
-  (7) F1 - G2, classes III - V, 2.60 < Hbeta < 2.72, 0.22 < (b-y)o < 0.39
-  (8) G2 - M2, classes  IV - V, 0.20 < m0    < 0.76, 0.39 < (b-y)o < 1.00
+  1. B0 - A0, classes III - V, 2.59 < Hbeta < 2.88,-0.20 <   c0   < 1.00
+  2. B0 - A0, class   Ia     , 2.52 < Hbeta < 2.59,-0.15 <   c0   < 0.40
+  3. B0 - A0, class   Ib     , 2.56 < Hbeta < 2.61,-0.10 <   c0   < 0.50
+  4. B0 - A0, class   II     , 2.58 < Hbeta < 2.63,-0.10 <   c0   < 0.10
+  5. A0 - A3, classes III - V, 2.87 < Hbeta < 2.93,-0.01 < (b-y)o < 0.06
+  6. A3 - F0, classes III - V, 2.72 < Hbeta < 2.88, 0.05 < (b-y)o < 0.22
+  7. F1 - G2, classes III - V, 2.60 < Hbeta < 2.72, 0.22 < (b-y)o < 0.39
+  8. G2 - M2, classes  IV - V, 0.20 < m0    < 0.76, 0.39 < (b-y)o < 1.00
 * `hbeta` (optional): H-beta line strength index. If it is not supplied, then by
   default its value will be `NaN` and the code will estimate a value based on by,
   m1,and c1. It is not used for stars in group 8.
@@ -280,20 +280,12 @@ Derive dereddened colors, metallicity, and Teff from Stromgren colors.
 
 ### Example ###
 
-Suppose 5 stars have the following Stromgren parameters
-
-by = [-0.001 ,0.403, 0.244, 0.216, 0.394]
-m1 = [0.105, -0.074, -0.053, 0.167, 0.186]
-c1 = [0.647, 0.215, 0.051, 0.785, 0.362]
-hbeta = [2.75, 2.552, 2.568, 2.743, 0]
-nn = [1,2,3,7,8]
-
-Determine the stellar parameters
+Determine the stellar parameters of 5 stars given their Stromgren parameters
 
 ```jldoctest
 julia> using AstroLib
 
-julia> by = [-0.001 ,0.403, 0.244, 0.216, 0.394];
+julia> by = [-0.001, 0.403, 0.244, 0.216, 0.394];
 
 julia> m1 = [0.105, -0.074, -0.053, 0.167, 0.186];
 
@@ -301,10 +293,10 @@ julia> c1 = [0.647, 0.215, 0.051, 0.785, 0.362];
 
 julia> hbeta = [2.75, 2.552, 2.568, 2.743, 0];
 
-julia> nn = [1,2,3,7,8];
+julia> nn = [1, 2, 3, 7, 8];
 
 julia> uvbybeta.(by, m1, c1, nn, hbeta)
-5-element Array{NTuple{5,Float64},1}:
+5-element Vector{NTuple{5, Float64}}:
  (13057.535222326893, -0.27375469585031265, 0.04954396423248884, -0.008292894218734928, 2.7136529525371897)
  (14025.053834219656, -6.907050783073221, 0.4140562248995983, NaN, 73.50771722263974)
  (18423.76405400214, -5.935816553877892, 0.2828247876690783, NaN, 39.84106215808709)
diff --git a/src/ydn2md.jl b/src/ydn2md.jl
index 6386dd3..647842c 100644
--- a/src/ydn2md.jl
+++ b/src/ydn2md.jl
@@ -30,7 +30,7 @@ Find the date of the 60th and 234th days of the year 2016.
 julia> using AstroLib
 
 julia> ydn2md.(2016, [60, 234])
-2-element Array{Dates.Date,1}:
+2-element Vector{Dates.Date}:
  2016-02-29
  2016-08-21
 ```
@@ -39,6 +39,5 @@ julia> ydn2md.(2016, [60, 234])
 
 `ymd2dn` converts from a date to day of the year.
 """
-function ydn2md(year::Integer, day::Integer)
-    return Dates.firstdayofyear(Date(year)) + Dates.Day(day - 1)
-end
+ydn2md(year::Integer, day::Integer) =
+    Dates.firstdayofyear(Date(year)) + Dates.Day(day - 1)
diff --git a/src/ymd2dn.jl b/src/ymd2dn.jl
index 42645de..01dd7a3 100644
--- a/src/ymd2dn.jl
+++ b/src/ymd2dn.jl
@@ -30,7 +30,7 @@ year).
 julia> using AstroLib, Dates
 
 julia> ymd2dn.([Date(2015, 3, 5), Date(2016, 3, 5)])
-2-element Array{Int64,1}:
+2-element Vector{Int64}:
  64
  65
 ```