Collection of Common Astronomical Conversion Routines and Functions

Description

Various functions for converting between commonly used coordinate systems in astronomy and making cosmological calculations.

Details

There are a number of functions included, but the most useful for astronomy conversions are the decimal degrees to DMS/HMS formats used at many telescopes: deg2dms, deg2hms, dms2deg, hms2deg. It also contains functions for various cosmological calculations (i.e. distance, volume and age for different cosmologies and redshifts).

Author(s)

Aaron Robotham

Maintainer: Aaron Robotham <aaron.robotham@uwa.edu.au>

Cosmology parameter data included in celestial package.

Description

cosref: Cosmology H0 / OmegaM / OmegaL / OmegaR (via OmegaM/zeq) and Sigma8 parameters taken from Planck (13/15/18), WMAP (1/3/5/9), Millennium Simulation and GiggleZ. Not all of these exist for each source, so NA values are used in these cases.

For Planck we use the second column of the main cosmology table, which does no use external data.

Usage

data(cosref)

Details

The included data is a table of the following cosmological parameters:

Author(s)

Aaron Robotham

References

See Also

cosvol, cosmap, cosdist, cosgrow

Examples

data(cosref)
cosref[cosref[,'Ref']=='Planck',]

IAU name creator.

Description

Creates IAU legal names for objects given coordinates, name and epoch.

Usage

IAUID(ra, dec, name = "GAMA", epoch = "J")

Arguments

Value

Text string that outputs an IAU legal name for an object.

Author(s)

Aaron Robotham

Examples

IAUID(123.45,67.89,'GAMA','J')

Sky matching

Description

These functions allows the user to match a reference set of sky coordinates against a comparison set of sky coordinates. The match radius can be varied per source (all matches per source are given within this radius), and mutual best matches are also extracted. coordmatch should be used for finding multiple matches and coordmatchsing should be used when trying to find matches around a single source. internalclean is a utility function that will remove closely duplicated objects via some tiebreak criterion, and is probably only of interest to advanced users trying to clean catalogues that were produced from overlapping frames.

Usage

coordmatch(coordref, coordcompare, rad = 2, inunitref = "deg", inunitcompare = "deg",
radunit = "asec", sep = ":", kstart = 10, ignoreexact = FALSE, ignoreinternal=FALSE,
matchextra = FALSE, smallapprox=FALSE)
coordmatchsing(RAref,Decref, coordcompare, rad=2, inunitref = "deg",
inunitcompare="deg", radunit='asec', sep = ":", ignoreexact=FALSE, smallapprox=FALSE)
internalclean(RA, Dec, rad=2, tiebreak, decreasing = FALSE, inunit="deg", radunit='asec',
sep = ":")

Arguments

Details

For coordmatch the main matching is done using nn2 that comes as part of the RANN package. coordmatch adds a large amount of sky coordinate oriented functionality beyond the simple implementation of nn2. For single object matches coordmatchsing should be used since it is substantially faster in this regime (making use of direct dot products).

ignoreexact is more strict in a sense since all objects exactly matching are ignored, whereas with ignoreinternal only identical row IDs are interpretted as being the same object.

Value

The output of coordmatch is a list containing:

The output of coordmatchsing is a list containing:

Author(s)

Aaron Robotham

See Also

hms2deg, dms2deg, sph2car

Examples

set.seed(666)

#Here we make objects in a virtual 1 square degree region

mocksky=cbind(runif(1e3), runif(1e3))

#Now we match to find all objects within an arc minute, ignoring self matches

mockmatches=coordmatch(mocksky, mocksky, ignoreexact=TRUE, rad=1, radunit='amin')

#Now we match to find all objects with varying match radii, ignoring self matches

mockmatchesvary=coordmatch(mocksky, mocksky, ignoreexact=TRUE, rad=seq(0,1,length=1e3),
radunit='amin')

#We can do this also by using the internal table match mode:

mockmatchesvary2=coordmatch(mocksky, rad=seq(0,1,length=1e3), radunit='amin')

#Check that this looks the same (should be identical with all zeroes):

summary(mockmatchesvary$bestmatch-mockmatchesvary2$bestmatch)

Transforms 3D cartesian coordinates to spherical coordinates

Description

Transforms 3D cartesian coordinates to spherical coordinates. The user can choose to return the spherical coordinates in degrees or radians.

Usage

car2sph(x, y, z, deg = TRUE)

Arguments

Details

This is a low level function that is used for plot transformations.

Value

A data.frame is returned containing the columns long (longitude), lat (latitude) and radius.

Author(s)

Aaron Robotham

See Also

sph2car

Examples

print(car2sph(x=1,y=1,z=0,deg=TRUE))

Navarro Frenk and White profile

Description

Density and total mass values for Navaro Frenk and White (NFW) profiles

Usage

cosNFW(Rad=0, Rho0=2.412e15, Rs=0.03253)
cosNFWmass_c(Rho0=2.412e15, Rs=0.03253, c=5, Munit = 1, Lunit = 1e+06)
cosNFWmass_Rmax(Rho0=2.412e15, Rs=0.03253, Rmax=0.16265, Munit = 1, Lunit = 1e+06)
cosNFWvcirc(Rad = 0.16264, Mvir = 1e+12, c = 5, f = Inf, z = 0, H0 = 100, OmegaM = 0.3,
OmegaL = 1 - OmegaM - OmegaR, OmegaR = 0, Rho = "crit", Dist = "Co", DeltaVir = 200,
Munit = 1, Lunit = 1e+06, Vunit = 1000, ref)
cosNFWvesc(Rad = 0.16264, Mvir = 1e+12, c = 5, f = Inf, z = 0, H0 = 100, OmegaM = 0.3,
OmegaL = 1 - OmegaM - OmegaR, OmegaR = 0, Rho = "crit", Dist = "Co", DeltaVir = 200,
Munit = 1, Lunit = 1e+06, Vunit = 1000, ref)
cosNFWsigma(Rad=0.03253, Rs=0.03253, c=5, z = 0, H0 = 100, OmegaM = 0.3,
OmegaL = 1-OmegaM-OmegaR, OmegaR=0, Rho = "crit", DeltaVir = 200, Munit = 1,
Lunit = 1e+06, Vunit = 1000, ref)
cosNFWsigma_mean(Rad=0.03253, Rs=0.03253, c=5, z = 0, H0 = 100, OmegaM = 0.3,
OmegaL = 1-OmegaM-OmegaR, OmegaR=0, Rho = "crit", DeltaVir = 200, Munit = 1,
Lunit = 1e+06, Vunit = 1000, ref)
cosNFWgamma(Rad=0.03253, Rs=0.03253, c=5, SigmaC=1, z = 0, H0 = 100,
OmegaM = 0.3, OmegaL = 1-OmegaM-OmegaR, OmegaR=0, Rho = "crit", DeltaVir = 200,
Munit = 1, Lunit = 1e+06, Vunit = 1000, ref)
cosNFWduffym2c(M=2e12, z = 0, H0 = 100, OmegaM = 0.3, OmegaL = 1-OmegaM-OmegaR,
OmegaR=0, Rho = "crit", A=6.71, B=-0.091, C=-0.44, Munit = 1, ref)

Arguments

Details

These functions calculate various aspects of the NFW profile.

Value

Author(s)

Aaron Robotham

References

Duffy A.R., et al., 2008, MNRAS, 390L

Navarro J.F., Frenk C.S., White Simon D.M., 1996, ApJ, 462

Wright C.O. & Brainerd T.G., 2000, ApJ, 534

See Also

cosvol, cosmap, cosdist, cosgrow, coshalo

Examples

#What difference do we see if we use the rad_mean200 radius rather than rad_crit200

rad_crit200=coshaloMvirToRvir(1e12,Lunit=1e6)
rad_mean200=coshaloMvirToRvir(1e12,Lunit=1e6,Rho='mean')
cosNFWmass_Rmax(Rmax=rad_crit200) #By construction we should get ~10^12 Msun/h
cosNFWmass_Rmax(Rmax=rad_mean200) #For the same profile this is a factor 1.31 larger

#Shear checks:

plot(10^seq(-2,2,by=0.1), cosNFWgamma(10^seq(-2,2,by=0.1),Rs=0.2,c=10), type='l',
log='xy', xlab='R/Rs', ylab='ESD')
legend('topright', legend=c('Rs=0.2','c=10'))

#How do critical, mean 200 and 500 masses evolve with redshift? Let's see:

zseq=10^seq(-2, 1, by=0.1)
con=seq(2, 20, by=0.01)
concol=rainbow(length(con), start=0, end=5/6)
rad_crit200=coshaloMvirToRvir(1, z=zseq, Rho='crit', DeltaVir=200, ref='Planck15')
rad_crit500=coshaloMvirToRvir(1, z=zseq, Rho='crit', DeltaVir=500, ref='Planck15')
rad_mean200=coshaloMvirToRvir(1, z=zseq, Rho='mean', DeltaVir=200, ref='Planck15')
rad_mean500=coshaloMvirToRvir(1, z=zseq, Rho='mean', DeltaVir=500, ref='Planck15')
rad_vir=coshaloMvirToRvir(1, z=zseq, Rho='crit', DeltaVir='get', ref='Planck15')

plot(1, 1, type='n', xlim=c(0.01,10), ylim=c(0.8,1.55), xlab='Redshift',
ylab='M200c / M500c', log='x')
for(i in 1:length(con)){
lines(zseq, cosNFWmass_Rmax(Rho0=1, Rs=rad_crit200[1]/con[i], Rmax=rad_crit200)/
cosNFWmass_Rmax(Rho0=1, Rs=rad_crit200[1]/con[i], Rmax=rad_crit500), col=concol[i])
}

plot(1, 1, type='n', xlim=c(0.01,10), ylim=c(0.8,1.55), xlab='Redshift',
ylab='M200m / M500m', log='x')
for(i in 1:length(con)){
lines(zseq, cosNFWmass_Rmax(Rho0=1, Rs=rad_crit200[1]/con[i], Rmax=rad_mean200)/
cosNFWmass_Rmax(Rho0=1, Rs=rad_crit200[1]/con[i], Rmax=rad_mean500), col=concol[i])
}

plot(1, 1, type='n', xlim=c(0.01,10), ylim=c(0.8,1.55), xlab='Redshift',
ylab='M200m / M200c',log='x')
for(i in 1:length(con)){
lines(zseq, cosNFWmass_Rmax(Rho0=1, Rs=rad_crit200[1]/con[i], Rmax=rad_mean200)/
cosNFWmass_Rmax(Rho0=1, Rs=rad_crit200[1]/con[i], Rmax=rad_crit200), col=concol[i])
}

plot(1, 1, type='n', xlim=c(0.01,10), ylim=c(0.8,1.55), xlab='Redshift',
ylab='M500m / M500c', log='x')
for(i in 1:length(con)){
lines(zseq, cosNFWmass_Rmax(Rho0=1, Rs=rad_crit200[1]/con[i], Rmax=rad_mean500)/
cosNFWmass_Rmax(Rho0=1, Rs=rad_crit200[1]/con[i], Rmax=rad_crit500), col=concol[i])
}

plot(1, 1, type='n', xlim=c(0.01,10), ylim=c(0.8,1.55), xlab='Redshift',
ylab='Mvir / M200c',log='x')
for(i in 1:length(con)){
lines(zseq, cosNFWmass_Rmax(Rho0=1, Rs=rad_crit200[1]/con[i], Rmax=rad_vir)/
cosNFWmass_Rmax(Rho0=1, Rs=rad_crit200[1]/con[i], Rmax=rad_crit200), col=concol[i])
}

plot(1, 1, type='n', xlim=c(0.01,10), ylim=c(0.8,1.55), xlab='Redshift',
ylab='Mvir / M200m',log='x')
for(i in 1:length(con)){
lines(zseq, cosNFWmass_Rmax(Rho0=1, Rs=rad_crit200[1]/con[i], Rmax=rad_vir)/
cosNFWmass_Rmax(Rho0=1, Rs=rad_crit200[1]/con[i], Rmax=rad_mean200), col=concol[i])
}

plot(zseq, rad_crit200/rad_crit500, type='l', xlim=c(0.01,10), ylim=c(0.8,1.55),
xlab='Redshift', ylab='R200 / R500', log='x')

plot(zseq, rad_mean200/rad_crit200, type='l', xlim=c(0.01,10), ylim=c(0.8,1.55), 
xlab='Redshift', ylab='Rm / Rc', log='x')

plot(zseq, rad_vir/rad_crit200, type='l', xlim=c(0.01,10), ylim=c(0.8,1.55), 
xlab='Redshift', ylab='Rvir / R200c', log='x')

plot(zseq, rad_vir/rad_mean200, type='l', xlim=c(0.01,10), ylim=c(0.8,1.55), 
xlab='Redshift', ylab='Rvir / R200c', log='x')

#R200m and R200c go either side of Rvir, so by cosmic conspiracy the mean is nearly flat:

plot(zseq, 2*rad_vir/(rad_mean200+rad_crit200), type='l', xlim=c(0.01,10),
ylim=c(0.8,1.55), xlab='Redshift', ylab='2Rvir / (R200c+R200m)', log='x')

#To check Vcirc and Vesc for a 10^12 Msun halo:

plot(0:400, cosNFWvcirc(0:400,f=1,Lunit=1e3), type='l', lty=1, xlab='R / kpc',
ylab='V / km/s', ylim=c(0,500))
lines(0:400, cosNFWvesc(0:400,f=1,Lunit=1e3), lty=2)
legend('topright', legend=c('Vel-Circ','Vel-Escape'), lty=c(1,2))
abline(v=coshaloMvirToRvir(Lunit=1e3), lty=3)

Cosmological distance calculations

Description

These functions allow comoving, angular size and luminosity distances to be calculated for a given redshift, it can also return look back time. They use curvature correctly, calculated internally using the relation OmegaM+OmegaL+OmegaR+OmegaK=1, but by default they assume a flat Universe where only OmegaM needs to be specified and OmegaR=0 (so no radiation pressure at any epoch).

Usage

cosdist(z=1, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0 = -1, wprime = 0,
age=FALSE, ref, error=FALSE)
cosdistz(z=1)
cosdistzeff(zref = 1, zem = 2)
cosdista(z=1)
cosdistCoDist(z=1, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0 = -1,
wprime = 0, ref)
cosdistLumDist(z=1, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0 = -1,
wprime = 0, ref)
cosdistAngDist(z=1, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0 = -1,
wprime = 0, ref)
cosdistAngDist12(z1=1, z2=2, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0,
w0 = -1, wprime = 0, ref)
cosdistCoDistTran(z=1, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0 = -1,
wprime = 0, ref)
cosdistCoDist12ang(z1=1, z2=2, ang=0, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR,
OmegaR=0, w0 = -1, wprime = 0, inunit='deg', ref) 
cosdistLumDist12ang(z1=1, z2=2, ang=0, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR,
OmegaR=0, w0 = -1, wprime = 0, inunit='deg', ref)
cosdistAngDist12ang(z1=1, z2=2, ang=0, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR,
OmegaR=0, w0 = -1, wprime = 0, inunit='deg', ref)
cosdistzem12ang(z1=1, z2=2, ang=0, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR,
OmegaR=0, w0 = -1, wprime = 0, inunit='deg', ref)
cosdistzeff12ang(z1=1, z2=2, ang=0, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR,
OmegaR=0, w0 = -1, wprime = 0, inunit='deg', ref)
cosdistDistMod(z=1, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0 = -1,
wprime = 0, ref)
cosdistAngScale(z=1, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0 = -1,
wprime = 0, ref)
cosdistAngSize(z=1, Size=1, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0=-1,
wprime=0, Dim=1, Dist='Co',  outunit='deg', ref)
cosdistAngArea(z=1, Size=1, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0=-1,
wprime=0, Dim=2, Dist='Co',  outunit='deg2', ref)
cosdistCoVol(z=1, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0 = -1,
wprime = 0, ref)
cosdistHubTime(H0=100)
cosdistUniAgeNow(z=1, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0 = -1,
wprime = 0, ref)
cosdistUniAgeAtz(z=1, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0 = -1,
wprime = 0, ref)
cosdistTravelTime(z=1, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0 = -1,
wprime = 0, ref)
cosdistRelError(z=1, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0 = -1, wprime = 0,
ref)
cosdistCrit(z_lens=1, z_source=2, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0,
w0 = -1, wprime = 0, ref)

Arguments

Details

Functions are largely based on D. W. Hogg et al. 1999 and Wright et al. 2006.

Negative value of z> -1 are allowed, which produces future predictions based on present day cosmology.

cosdistAngDist12 is only available for OmegaK>=0.

Value

If age=TRUE is set then additional age-related information is calculated for each z as extra columns:

If error=TRUE is set then the relative error for distance/age values is calculated for each z as an extra column:

The outputs of the standalone functions are:

Author(s)

Aaron Robotham

References

Based on the equations in:

Davis T.M. & Lineweaver, Charles H., 2004, PASA, 21, 97

Hogg D.W., 1999, arXiv, 9905116

Liske J., 2000, MNRAS, 319, 557L

Peacock J.A., 1999, Cosmological Physics, Cambridge University Press

Wright E.L., 2006, PASP, 118, 1711

See Also

cosvol, cosmap, cosgrow, cosref, cosNFW

Examples

## Not run: 
cosdist(0.3,70,age=TRUE)
cosdist(0.3,70,age=TRUE,ref='Planck')
cosdistz(0.3)
cosdista(0.3)
cosdistCoDist(0.3,70)
cosdistLumDist(0.3,70)
cosdistAngDist(0.3,70)
cosdistAngDist12(0.3,0.5,70)
cosdistCoDistTran(0.3,70)
cosdistCoDist12ang(0,2,10)
cosdistDistMod(0.3,70)
cosdistAngScale(0.3,70)
cosdistAngSize(0.3,1,70)
cosdistCoVol(0.3,70)
cosdistHubTime(70)
cosdistUniAgeNow(0.3,70)
cosdistUniAgeAtz(0.3,70)
cosdistTravelTime(0.3,70)
cosdistRelError(0.3)
cosdistCrit(0.3,0.5,70)
cosdistzeff(1,2)
cosdistzem12ang(1,2)
cosdistzeff12ang(1,2)

#A check of the comoving separation between objects function:

cosdistCoDistTran(2,OmegaM = 0.3, OmegaL=1)*sin(pi/180)
cosdistCoDist12ang(2,2,ang=1,OmegaM=0.3,OmegaL=1)

#Very close, however cosdistCoDist12ang lets us go further:

cosdistCoDist12ang(1,2,ang=10,OmegaM=0.3,OmegaL=1)
cosdistCoDist12ang(2,2,ang=180,OmegaM=0.3,OmegaL=1)

#The second number should be be the same as:

cosdistCoDist(2,OmegaM=0.3,OmegaL=1)*2

#Example 1 by John Peacock for EDS Universe (answer should be nearly 3):

cosdistzem12ang(3,4,56.4,H0=100,OmegaM=1,OmegaL=0)

#Example 2 by John Peacock for EDS Universe (answer should be nearly 2995 Mpc/h):

cosdistCoDist12ang(3,4,56.4,H0=100,OmegaM=1,OmegaL=0)

#Example 3 by John Peacock for Milne Universe (answer should be nearly 5294 Mpc/h):

cosdistCoDist12ang(3,4,56,H0=100,OmegaM=0,OmegaL=0)

#Example 4 by John Peacock for Milne Universe (answer should be nearly 4.846):

cosdistzeff12ang(3,4,56,H0=100,OmegaM=0,OmegaL=0)

#Example 5 by John Peacock for Milne Universe (answer should be nearly 364 Mpc/h):

cosdistAngDist12ang(3,4,56,H0=100,OmegaM=0,OmegaL=0)

#Nice plot of distance estimates:

redshifts=seq(0,3,by=0.01)
plot(redshifts, cosdistCoDist(redshifts, ref='planck'), type='l', col='darkgreen',
xlab='Redshift / z', ylab='Distance / Mpc')
lines(redshifts, cosdistLumDist(redshifts, ref='planck'), col='red')
lines(redshifts, cosdistAngDist(redshifts, ref='planck'), col='blue')
legend('topleft', legend=c('Comoving Distance', 'Luminosity Distance', 'Angular Diameter Distance'),
col=c('darkgreen', 'red', 'blue'),lty=1)

plot(redshifts, cosdistTravelTime(redshifts, ref='planck'), type='l',
xlab='Redshift / z', ylab='Light travel time / Yrs')

#Actual time example (Figure 1 of Davis & Lineweaver 2004)
zseq=10^seq(-2,6,len=1e3)-1
dists=cosdistCoDist(zseq, ref='737')*0.00326
times=cosdistTravelTime(zseq, ref='737')
plot(dists, times, type='l', xlab='Comoving Distance / Glyr',
ylab='Time / Gyr')
abline(v=0, h=0, lty=1)
abline(h=c(min(times), max(times)), lty=2)
abline(v=c(min(dists), max(dists)), lty=2)

#Conformal time example (Figure 1 of Davis & Lineweaver 2004):
#Mpc to Glyr conversion is 0.00326

zseq=10^seq(-2,6,len=1e3)-1
dists=cosdistCoDist(zseq, ref='737')*0.00326
plot(dists, dists, type='l',
xlab='Comoving Distance / Glyr', ylab='Conformal Time / Gyr')
abline(v=0, h=0, lty=1)
abline(h=c(min(dists), max(dists)), lty=2)
abline(v=c(min(dists), max(dists)), lty=2)

## End(Not run)

Cosmological growth and evolution calculations

Description

These functions allow various properties of the expansion of the Universe to be calculated: e.g. OmegaM/OmegaL/OmegaR/OmegaK for ay redshift, growth rate and growth factor, sigma8, and RhoCrit and RhoMean. They use curvature correctly, calculated internally using the relation OmegaM+OmegaL+OmegaR+OmegaK=1, but by default they assume a flat Universe where only OmegaM needs to be specified and OmegaR=0 (so no radiation pressure at any epoch).

Usage

  cosgrow(z=1, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0 = -1, wprime = 0,
  Sigma8=0.8, fSigma8=FALSE, Dist='Co',
  Mass='Msun', ref)
  cosgrowz(z = 1)
  cosgrowa(z = 1)
  cosgrowH(z=1, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0=-1, wprime=0,
  ref)
  cosgrowCoVel(z=1, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0=-1,
  wprime=0, ref)
  cosgrowPecVel(z=1, zob=1)
  cosgrowOmegaM(z=1, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0=-1, wprime=0, ref)
  cosgrowOmegaL(z=1, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0=-1, wprime=0, ref)
  cosgrowOmegaR(z=1, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0=-1, wprime=0, ref)
  cosgrowOmegaK(z=1, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0=-1, wprime=0, ref)
  cosgrowDecelq(z=1, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0=-1, wprime=0, ref)
  cosgrowEoSwDE(z=1, w0=-1, wprime=0)
  cosgrowRhoDE(z=1,w0=-1, wprime=0, rhoDE=1)
  cosgrowFactor(z=1, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0=-1, wprime=0, ref)
  cosgrowRate(z=1, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0=-1, wprime=0,
  Sigma8=0.8, fSigma8=FALSE, ref)
  cosgrowSigma8(z=1, OmegaM=0.3,OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0=-1, wprime=0,
  Sigma8=0.8, ref)
  cosgrowFactorApprox(z=1, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0=-1, wprime=0,
  ref)
  cosgrowRateApprox(z=1, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0=-1, wprime=0,
  Sigma8=0.8, fSigma8=FALSE, ref)
  cosgrowSigma8Approx(z=1, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0=-1, wprime=0,
  Sigma8=0.8, ref)
  cosgrowRhoCrit(z=1, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0=-1,
  wprime=0, Dist='Co', Mass='Msun', ref)
  cosgrowRhoMean(z=1, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, w0=-1,
  wprime=0, Dist='Co', Mass='Msun', ref)
  cosgrowDeltaVir(z=1, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0, ref)

Arguments

Details

The above functions are heavily based on the equations in Hamilton A.J.S., 2001, MNRAS 322 419 and Lahav O., et al., 1991, MNRAS, 251, 136.

Negative value of z> -1 are allowed, which produces future predictions based on present day cosmology.

The approximation routines are generally accurate to sub 1 percent, and since they do not involve numerical integration they are substantially faster when computing large grids of numbers, i.e. they are recommended for plots, since the accuracy is sub the line width.

Value

The outputs of the standalone functions are:

Note

The difference between RhoCrit and RhoMean at z=0 is simply RhoMean=RhoCrit*OmegaM. Corrected for 1/(1+z)^3 RhoMean stays constant with redshift (as it should- to first order we do not gain or lose mass within a comoving volume).

The growth rate and growth factor does not make use of OmegaR in the cosgrow function, hence OmegaR cannot be provided in the individual functions. This is because correctly accounting for the effect of radiation pressure before the surface of last scattering (z~1100) on the growth rate of structure is highly complex, and beyond the scope of this package. In the case of cosgrow, even if OmegaR is specified it is, in effect, set to zero when making growth factor and rate calculations.

The evolution of the dark matter equation of state (w) is parameterised as described in Wright (2006).

It is important to remember that H is in physical units for both the numerator and denominator (i.e. 'proper' at a given redshift, so the units are km/s / pMpc). To ask the question "is the Universe accelerating?" is to really ask "is the expansion factor accelerating?". This requires the denominator to be in comoving units (so rescaling for proper distances today) and measuring the differential with time or redshift. You will only find an accelerating Universe when dividing H(z)/(1+z)! See the examples to see how we can find this location, and that it is consistent with the start of acceleration calculated from the decleration parameter (q) directly.

Author(s)

Aaron Robotham

References

Based on the equations in:

Bryan & Norman, 1998, ApJ, 495, 80

Davis T.M. & Lineweaver, Charles H., 2004, PASA, 21, 97

Davis T.M. & Scrimgeour M.I., 2014, MNRAS, 442, 1117

Hamilton A.J.S., 2001, MNRAS 322 419

Lahav O., et al., 1991, MNRAS, 251, 136

Peacock J.A., 1999, Cosmological Physics, Cambridge University Press

Wright E.L., 2006, PASP, 118, 1711

See Also

cosvol, cosmap, cosdist, cosref, coshalo

Examples

  cosgrow(0.3)
  cosgrow(0.3,ref='Planck')
  cosgrowz(0.3)
  cosgrowa(0.3)
  cosgrowH(0.3)
  cosgrowCoVel(0.3)
  cosgrowPecVel(0.3,0.31)
  cosgrowOmegaM(0.3)
  cosgrowOmegaL(0.3)
  cosgrowOmegaK(0.3)
  sum(cosgrowOmegaM(0.3)+cosgrowOmegaL(0.3)+cosgrowOmegaK(0.3)) #Still 1.
  cosgrowDecelq(0.3)
  cosgrowEoSwDE(0.3)
  cosgrowFactor(0.3)
  cosgrowFactorApprox(0.3) #Approximation better than 1% for reasonable cosmologies.
  cosgrowRate(0.3)
  cosgrowRateApprox(0.3) #Approximation better than 1% for reasonable cosmologies.
  cosgrowRhoCrit(0.3)
  cosgrowRhoMean(0.3)
  cosgrowRhoMean(0)-cosgrowRhoMean(2,Dist='Ang')/(1+2)^3 #Mass is conserved in co-vol
  cosgrowRhoMean(0)-cosgrowRhoMean(10,Dist='Co') #Mass is conserved in co-vol
  
  # Various recessional velocities (see Figure 2 of Davis & Lineweaver 2004):
  
  plot(10^seq(-1,4,by=0.01), cosgrowCoVel(10^seq(-1,4,by=0.01), ref='planck')
  /299792.458, type='l', log='x', xlab='z', ylab='Cosmological Recession Velocity / c')
  lines(10^seq(-1,4,by=0.01), cosgrowPecVel(0,10^seq(-1,4,by=0.01))/299792.458, col='red')
  lines(10^seq(-1,4,by=0.01), 10^seq(-1,4,by=0.01), col='blue')
  abline(h=1,v=1.5,lty=2)
  legend('topleft', legend=c('GR', 'SR', 'Approx (cz)', 'Superluminal'), lty=c(1,1,1,2),
  col=c('black','red','blue','black'))
  
  # Comparison of the various energy densities that make up the Universe for Planck 2013:
  
  plot(cosdistUniAgeAtz(10^seq(-3,4.9,by=0.1), ref='Planck')*1e9,
  cosgrowRhoCrit(10^seq(-3,4.9,by=0.1), ref='Planck', Dist='m', Mass='kg')*
  cosgrowOmegaR(10^seq(-3,4.9,by=0.1), ref='Planck'), type='l',log='xy',
  xlab='Years since Universe formed', ylab=expression('Energy Density'*(kg/m^3)))
  
  lines(cosdistUniAgeAtz(10^seq(-3,4.9,by=0.1), ref='Planck')*1e9,
  cosgrowRhoCrit(10^seq(-3,4.9,by=0.1), ref='Planck', Dist='m', Mass='kg')*
  cosgrowOmegaM(10^seq(-3,4.9,by=0.1), ref='Planck'), col='red')
  
  lines(cosdistUniAgeAtz(10^seq(-3,4.9,by=0.1), ref='Planck')*1e9,
  cosgrowRhoCrit(10^seq(-3,4.9,by=0.1), ref='Planck', Dist='m', Mass='kg')*
  cosgrowOmegaL(10^seq(-3,4.9,by=0.1), ref='Planck'), col='blue')

  abline(v=cosdistUniAgeAtz(0.33,ref='Planck')*1e9,lty=2) # Matter = Vacuum
  abline(v=cosdistUniAgeAtz(3391,ref='Planck')*1e9,lty=2) # Matter = Radiation
  
  legend('topright', legend=c('Radiation Energy Density', 'Matter Energy Density',
  'Vacuum Energy Density'), lty=1, col=c('black','red','blue'))
  
  # Where's the acceleration?
  plot(seq(0,2,len=1e3),cosgrowH(seq(0,2,len=1e3)),type='l',xlab='z',
  ylab='H(z) / km/s / pMpc')
  # There it is!
  plot(seq(0,2,len=1e3),cosgrowH(seq(0,2,len=1e3))/(1+seq(0,2,len=1e3)),
  type='l',xlab='z',ylab='H(z) / km/s / cMpc')
  #When does it start accelerating?
  accel.loc=which.min(abs(cosgrowDecelq(seq(0,2,len=1e3))))
  abline(v=seq(0,2,len=1e3)[accel.loc],lty=2)

Virial halo conversion functions

Description

All 6 Virial parameter conversion functions. Each can map precisely to the other as a one paramter function.

Usage

coshaloMvirToSigma(Mvir=1e+12, z=0, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR,
OmegaR=0, Rho='crit', Dist='Co', DeltaVir=200, Munit=1, Lunit=1e6, Vunit=1e3, Dim=3, ref)
coshaloSigmaToMvir(Sigma=230, z=0, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR,
OmegaR=0, Rho='crit', Dist='Co', DeltaVir=200, Munit=1, Lunit=1e6, Vunit=1e3, Dim=3, ref)
coshaloMvirToRvir(Mvir=1e12, z=0, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR,
OmegaR=0, Rho='crit', Dist='Co', DeltaVir=200, Munit=1, Lunit=1e6, Vunit=1e3, Dim=3, ref)
coshaloRvirToMvir(Rvir=162.635, z=0, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR,
OmegaR=0, Rho='crit', Dist='Co', DeltaVir=200, Munit=1, Lunit=1e6, Vunit=1e3, Dim=3, ref)
coshaloSigmaToRvir(Sigma=230, z=0, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR,
OmegaR=0, Rho='crit', Dist='Co', DeltaVir=200, Munit=1, Lunit=1e6, Vunit=1e3, Dim=3, ref)
coshaloRvirToSigma(Rvir=162.635, z=0, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR,
OmegaR=0, Rho='crit', Dist='Co', DeltaVir=200, Munit=1, Lunit=1e6, Vunit=1e3, Dim=3, ref)
coshaloSigmaToTvir(Sigma=230, Vunit=1e3, Tunit='K', Type='halo', Dim=3) 

Arguments

Details

These functions allow for various analytic conversions between the 3 major properties related to virial radius: the mass, velocity dispresion and size. The default properties calculate properties for 1e12 Msun halos and assume masses in Msun, velocities in km/s and distances in Kpc.

Value

Author(s)

Aaron Robotham, Chris Power

References

coshaloSigmaToTvir based on the equations in:

Bryan & Norman, 1998, ApJ, 495, 80

See Also

cosvol, cosmap, cosdist, cosgrow, cosNFW

Examples

coshaloMvirToSigma(1e13) # Velocity in km/s
coshaloMvirToSigma(1e13, Vunit=1) # Velocity in m/s
coshaloSigmaToMvir(coshaloMvirToSigma(1e13, Vunit=1),Vunit=1)
coshaloMvirToRvir(1e13) #Radius in kpc
coshaloSigmaToRvir(coshaloMvirToSigma(1e13, Vunit=1),Vunit=1)

#Some sanity checks

rho_crit200=cosgrowRhoCrit(z=0)*200 #200 times rho critical at z=0
rho_mean200=cosgrowRhoMean(z=0)*200 #200 times rho mean at z=0
#For a 10^12 Msun/h halo, the radius in Mpc/h where the contained density equals rho_crit*200
rad_crit200=(1e12/rho_crit200*3/4/pi)^(1/3)
coshaloMvirToRvir(1e12,Lunit=1e6)-rad_crit200 # ~0 as expected
#For a 10^12 Msun/h halo, the radius in Mpc/h where the contained density equals rho_crit*200
rad_mean200=(1e12/rho_mean200*3/4/pi)^(1/3) # ~0 as expected
coshaloMvirToRvir(1e12,Lunit=1e6,Rho='mean')-rad_mean200

Cosmological Mapping Functions

Description

Functions for mapping from one arbitrary cosmological parameter to another. This includes the provision of a generic interpolation function and another exact value lookup.

Usage

cosmapval(val=50, cosparam="CoVol", H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0,
w0=-1, wprime=0, Sigma8=0.8, fSigma8=FALSE, zrange=c(-0.99,100), res=100, iter=8,
out='cos', degen='lo', ref)
cosmapfunc(cosparamx="CoVol", cosparamy="z", H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR,
OmegaR=0, w0=-1, wprime=0, Sigma8=0.8, fSigma8=FALSE, zrange=c(0,20), step='z', res=100,
degen='lo', ref)

Arguments

Details

The default zrange and res should be sufficient for most reasonable cosmologies if the approximate redshift location of the region to be mapped is entirely unknown.

Predictions into the future are possible if val is set to negative (distance and volume parameters) or below their present day value (age and growth parameters). However, many potential values are outside of the asymptotic limits, e.g. using the default 737 cosmology H is tending to 83.666, i.e. it will fail if you request H=83 but work if your resuest H=84.

The default res and iter for cosmapval is appropriate for most mappings with -0.99 < z < 100 using a fiducial 737 cosmology. If this proves insufficient (this should be obvious from error column) then increase both of these. Overall accuracy goes as res^iter.

Value

If out='cos', cosmapval contains the concatenation of the cosdist (with age=TRUE and error=TRUE) and cosgrow functions for parameter 'cosparam' at value 'val'. The 'z' and 'a' columns are only included once (from the output of cosdist). See cosdist and cosgrow for information on the cosdist and cosgrow outputs. If out='z', then cosmapval merely returns the corresponding redshifts.

The cosmapval output (when out='cos') includes an additional final column named 'MapError' which gives the approximate relative error of the values returned compared to the desired lookup location. Smaller is obviously better, but at the cost of computational time.

cosmapfunc uses base R approxfun to map cosparamx onto cosparamy between zrange[1] and zrange[2] in uniform steps of expansion factor (a=1/(1+z)). cosmofunc returns the output function created by approxfun.

Author(s)

Aaron Robotham

References

Based on the equations in:

Hogg D.W., 1999, arXiv, 9905116

Wright E.L., 2006, PASP, 118, 1711

See Also

cosdist, cosvol, cosgrow

Examples

## Not run: 
tempfunc=cosmapfunc('CoVol', 'UniAgeAtz')
tempfunc(50)

cosmapval(50:60, 'CoVol')

#A future prediction:

cosmapval(59, 'H', H0=70)

## End(Not run)

Orbital functions

Description

A variety of obital analysis functions. These are useful for setting up initial conditions for merging systems etc.

Usage

  cosorbVisViva(M=1e12, Rad=162.635, SemiMajRad=162.635, Munit=1, Lunit=1e3, Vunit=1)
  cosorbFreeFall(M1=1e12, M2=1, Rad=162.635, Munit=1, Lunit=1e3, Vunit=1, Tunit=1e9)
  cosorbRocheRad(M1=1e12, M2=1e10, Size=35.03865, Rfac=2.44)
  cosorbRocheSize(M1=1e12, M2=1e10, Rad=396.8294, Rfac=2.44)

Arguments

Details

These functions allow for various analytic conversions between the 3 major properties related to virial radius: the mass, velocity dispresion and size. The default properties calculate properties for 1e12 Msun halos and assume masses in Msun, velocities in km/s and distances in Kpc.

Value

cosorbVisViva function gives the required velocity in units of Vunit to create the specified orbit.

cosorbFreeFall function gives the free fall time to static initial velocity separated bodies..

cosorbRocheRad function gives the orbital radius at which the secondary will become stripped within a specified bound radius.

cosorbRocheSize function gives the limiting bound radius of the secondary for a given system.

Author(s)

Aaron Robotham, Chris Power

See Also

cosvol, cosmap, cosdist, cosgrow

Examples

  cosorbVisViva(M=1e15, Rad=1, Lunit=1e6)
  cosorbFreeFall(M1=1e15, M2=1, Rad=1, Lunit=1e6)
  cosorbRocheRad(M1=1e12, M2=1e12, Size=162.635, Rfac=2.44)
  cosorbRocheSize(M1=1e12, M2=1e12, Rad=396.8294, Rfac=2.44)

Driver & Robotham (2010) cosmic variance calculator

Description

The main cosmic variance calculator function taken from Driver & Robotham (2010). cosvarcar is an interface to the Cartesian coordinate version, whilst cosvarsph is a utility interface to give approximate cosmic variance for astronomy survey regions (usually defined by RA, Dec and redshift limits).

Usage

cosvarcar(aside = 50, bside = 50, cside = 50, regions = 1)
cosvarsph(long = c(129, 141), lat = c(-2, 3), zmax = 1, zmin = 0, regions = 1,
inunit='deg', sep=":")
cosvararea(area=60, zmax=1, zmin=0, regions=1, inunit='deg2')

Arguments

Details

These functions use the empircally motivated cosmic variance percentage formula provided in Driver & Robotham (2010) Eqn 4.

cosvarsph is a 'best effort' approximation of the comoving box subtended by the specified spherical coordinates using the following conversions:

CoDistLow = cosdistCoDist(z=zmin,H0=70,OmegaM=0.3)

CoDistHigh = cosdistCoDist(z=zmax,H0=70,OmegaM=0.3)

cside=CoDistHigh-CoDistLow

area=skyarea(long = long, lat = lat, inunit = inunit, outunit='deg2')[1]

volume=cosvol(area=area, zmax = zmax, zmin=zmin, H0 = 70, OmegaM = 0.3, inunit='deg2')[1]

aside=cos(mean(lat)*pi/180)*(abs(diff(long))/360)*(CoDistLow+cside/2)

bside=(abs(diff(long))/180)*(CoDistLow+cside/2)

scale=sqrt(volume*1e9/(aside*bside*cside))

aside=aside*scale

bside=bside*scale

return(cosvarcar(aside=aside, bside=bside, cside=cside, subsets=subsets))

cosvararea is a simplifed version of cosvarsph, where the assumption is that aside=bside (so the aspect ratio on the sky is 1:1).

Value

The output is the approximate percentage cosmic (or sample) variance that is expected for the volume specified.

Note

Many people get upset at the term 'cosmic variance' and prefer 'sample variance'. Whilst I am sympathetic to the argument, more astronomers are familiar with the former term.

These cosmic variance estimates are defined using SDSS at z~0.1, caveats abound at higher redshifts, but these numbers should serve as a reasonably conservative (i.e. pessimistic) upper limit.

Author(s)

Aaron Robotham and Simon Driver

References

Driver S.P. & Robotham A.S.G., 2010, MNRAS, 407, 2131

See Also

cosvol, skyarea

Examples

#Approximate CV of the GAMA equatorial regions:
cosvarsph(long=12, lat=5, zmax=0.5)*1/sqrt(3)
#Or using the GAMA sexigesimal coordinates (should be the same):
cosvarsph(long = c('11:36:0','12:24:0'), lat = c('-2:0:0','3:0:0'), zmax=0.5,
inunit='sex')*1/sqrt(3)
#Approximate CV of the SDSS:
cosvarsph(long=150, lat=100, zmax=0.3)

Cosmological volume calculator

Description

Given the sky area, two redshifts and the cosmology, this function calculates the comoving volume.

Usage

cosvol(area=60, zmax=1, zmin=0, H0=100, OmegaM=0.3, OmegaL=1-OmegaM-OmegaR, OmegaR=0,
w0=-1, wprime=0, inunit = "deg2", ref)

Arguments

Value

A 3 element vector. The first element (voltot) specifies the comoving volume of the requested cone segment in Gpc^3, the second element (volmeanz) specifies the mean redshift when mass is uniformly distributed in the volume, the third element (volmedz) specifies the median redshift when mass is uniformly distributed in the volume.

Author(s)

Aaron Robotham

References

Based on the equations in:

Hogg D.W., 1999, arXiv, 9905116

Wright E.L., 2006, PASP, 118, 1711

See Also

cosdist,skyarea, cosmap, cosgrow

Examples

#Approximate volume of the GAMA survey (area given in skyarea example, zmax is approx
#limit of main galaxy sample):
TotalGAMAvol=cosvol(293.82,0.6)[1]
print(paste('The GAMA survey volume is ~',round(TotalGAMAvol,2),'Gpc^3'))

#Approximate volume of SDSS (area given for DR7, zmax is approx limit of main galaxy sample):
TotalSDSSvol=cosvol(8423,0.3)[1]
print(paste('The SDSS survey volume is ~',round(TotalSDSSvol,2),'Gpc^3'))

#Change of reference cosmology
cosvol(293.82,0.6,ref='Planck')

Convert decimal degrees to dms format.

Description

Convert decimal degrees to dms (degrees, minutes, seconds) format. This is probably most useful for declination conversion, since dms is fairly standard method of presenting declination coordinates. The decimal degrees=d+m/60+s/3600. Degrees should range from -90 to +90.

Usage

deg2dms(deg, type='mat', sep=':', digits=2)

Arguments

Value

A data.frame with the columns degrees, minutes and seconds if type='mat'. If type='cat' then a vector of strings with separators defined by the 'sep' argument.

Author(s)

Aaron Robotham

See Also

dms2deg

Examples

print(deg2dms(12.345))
print(deg2dms(12.345,type='cat',sep=':'))
print(deg2dms(12.345,type='cat',sep='dms'))
print(deg2dms(12.345,type='cat',sep='DMS'))

Convert decimal degrees to hms format.

Description

Convert decimal degrees to hms (hours, minutes, seconds) format. This is probably most useful for right-ascension (RA) conversion, since hms is fairly standard method of presenting RA coordinates. The decimal degrees=15*h+15*m/60+15*s/3600 (i.e. there are 24 hours in 360 degrees). Degrees should range from 0 to 360.

Usage

deg2hms(deg, type='mat', sep=':', digits=2)

Arguments

Value

A data.frame with the columns degrees, minutes and seconds if type='mat'. If type='cat' then a vector of strings with separators defined by the 'sep' argument.

Author(s)

Aaron Robotham

See Also

hms2deg

Examples

deg2hms(123.456)
deg2hms(123.456,type='cat',sep=':')
deg2hms(123.456,type='cat',sep='hms')
deg2hms(123.456,type='cat',sep='HMS')

Convert DMS to degrees format.

Description

Convert DMS (degrees, minutes, seconds) to degrees format. This is probably most useful for declination conversion, since dms is fairly standard method of presenting declination coordinates. The decimal degrees=d+m/60+s/3600. Degrees should range from -90 to +90. Degrees and minutes should be integer and seconds can be decimal.

Usage

dms2deg(d,m,s,sign='d',sep=':')

Arguments

Value

A value of decimal degrees.

Author(s)

Aaron Robotham

See Also

deg2dms

Examples

print(dms2deg(70,45,19,-1))
print(dms2deg('-70:45:19'))
print(dms2deg('-70d45m19s',sep='dms'))
print(dms2deg(c('-70D45M19S','3D5M15S'),sep='DMS'))

Get Pixel Scale

Description

Given a FITSio of astro header, calculate the image pixel scale.

Usage

getpixscale(header, CD1_1 = 1, CD1_2 = 0, CD2_1 = 0, CD2_2 = 1)

Arguments

Details

In most cases users will simply provide a valid header to find the WCS, but you can enter the CD values explicitly. Calculating the pixel scale from the latter is almost trivial, but the option is there for the curious/lazy.

Value

Numeric scalar; the image pixscale in asec/pixel (so typically a value of 0.1-0.5 for modern survey instruments).

Author(s)

Aaron Robotham

Examples

## Not run: 
#The answer should be almost exactly 0.2 asec/pixel:

#Using FITSio and ProFit packages
image = readFITS(system.file("extdata", 'KiDS/G266035fitim.fits', package="ProFit"))
getpixscale(image$hdr)
#Using astro package
image = read.fits(system.file("extdata", 'KiDS/G266035fitim.fits', package="ProFit"))
getpixscale(image$hdr[[1]])

## End(Not run)

Convert hms to degrees format.

Description

Convert hms (hours, minutes, seconds) to degrees format. This is probably most useful for right ascension (RA) conversion, since hms is fairly standard method of presenting RA coordinates. The decimal degrees=15*h+15*m/60+15*s/3600. Should range between 0 and 24 hours. Hours and minutes should be integer and seconds can be decimal.

Usage

hms2deg(h,m,s,sep=':')

Arguments

Value

A value of decimal degrees.

Author(s)

Aaron Robotham

See Also

deg2hms

Examples

hms2deg(12,10,36)
hms2deg('12:10:36')
hms2deg('12h10m36s',sep='hms')
hms2deg(c('12H10M36S','3H4M10S'),sep='HMS')

Planck's Law and Related Functions

Description

Functions related to Planck's Law of thermal radiation.

Usage

cosplanckLawRadFreq(nu,Temp=2.725)
cosplanckLawRadWave(lambda,Temp=2.725)
cosplanckLawEnFreq(nu,Temp=2.725)
cosplanckLawEnWave(lambda,Temp=2.725)
cosplanckLawRadFreqN(nu,Temp=2.725)
cosplanckLawRadWaveN(lambda,Temp=2.725)
cosplanckPeakFreq(Temp=2.725)
cosplanckPeakWave(Temp=2.725)
cosplanckSBLawRad(Temp=2.725)
cosplanckSBLawRad_sr(Temp=2.725)
cosplanckSBLawEn(Temp=2.725)
cosplanckLawRadPhotEnAv(Temp=2.725)
cosplanckLawRadPhotN(Temp=2.725)
cosplanckCMBTemp(z,Temp=2.725)

Arguments

Details

The functions with Rad in the name are related the spectral radiance form of Planck's Law (typically designated I or B), whilst those with En are related to the spectral energy density form of Planck's Law (u), where u=4\pi I/c.

To calculate the number of photons in a mode we simply use E=h\nu=h c / \lambda.

Below h is the Planck constant, k_B is the Boltzmann constant, c is the speed-of-light in a vacuum and \sigma is the Stefan-Boltzmann constant.

cosplanckLawRadFreq is the spectral radiance per unit frequency version of Planck's Law, defined as:

B_\nu(\nu,T) = I_\nu(\nu,T) = \frac{2 h \nu^3}{c^2} \frac{1}{e^{h \nu / k_B T}-1}

cosplanckLawRadWave is the spectral radiance per unit wavelength version of Planck's Law, defined as:

B_\lambda(\lambda,T) = I_\lambda(\lambda,T) = \frac{2 h c^2}{\lambda^5} \frac{1}{e^{h c / \lambda k_B T}-1}

cosplanckLawRadFreqN is the number of photons per unit frequency, defined as:

B_\nu(\nu,T) = I_\nu(\nu,T) = \frac{2 \nu^2}{c^2} \frac{1}{e^{h \nu / k_B T}-1}

cosplanckLawRadWaveN is the number of photons per unit wavelength, defined as:

B_\lambda(\lambda,T) = I_\lambda(\lambda,T) = \frac{2 c}{\lambda^4} \frac{1}{e^{h c / \lambda k_B T}-1}

cosplanckLawEnFreq is the spectral energy density per unit frequency version of Planck's Law, defined as:

u_\nu(\nu,T) = \frac{8 \pi h \nu^3}{c^3} \frac{1}{e^{h\nu/k_B T}-1}

cosplanckLawEnWave is the spectral energy density per unit wavelength version of Planck's Law, defined as:

u_\lambda(\lambda,T) = \frac{8 \pi h c}{\lambda^5} \frac{1}{e^{h c / \lambda k_B T}-1}

cosplanckPeakFreq gives the location in frequency of the peak of I_\nu(\nu,T), defined as:

\nu_{peak} = 2.821 k_B T

cosplanckPeakWave gives the location in wavelength of the peak of I_\lambda(\lambda,T), defined as:

\lambda_{peak} = 4.965 k_B T

cosplanckSBLawRad gives the emissive power (or radiant exitance) version of the Stefan-Boltzmann Law, defined as:

j^* = \sigma T^4

cosplanckSBLawRad_sr gives the spectral radiance version of the Stefan-Boltzmann Law, defined as:

L = \sigma T^4/\pi

cosplanckSBLawEn gives the energy density version of the Stefan-Boltzmann Law, defined as:

\epsilon = 4 \sigma T^4 / c

Notice that j^* and L merely differ by a factor of \pi, i.e. L is per steradian.

cosplanckLawRadPhotEnAv gives the average energy of the emitted black body photon, defined as:

<E_{phot}> = 3.729282 \times 10^{-23} T

cosplanckLawRadPhotN gives the total number of photons produced by black body per metre squared per second per steradian, defined as:

N_{phot} = 1.5205 \times 10^{15} T^3 / \pi

Various confidence building sanity checks of how to use these functions are given in the Examples below.

Value

Planck's Law in terms of spectral radiance:

Planck's Law in terms of spectral energy density:

Photon counts:

Peak locations (via Wien's displacement law):

Stefan-Boltzmann Law:

Photon properties:

Cosmic Microwave Background:

Author(s)

Aaron Robotham

References

Marr J.M., Wilkin F.P., 2012, AmJPh, 80, 399

See Also

cosgrow

Examples

#Classic example for different temperature stars:

waveseq=10^seq(-7,-5,by=0.01)
plot(waveseq, cosplanckLawRadWave(waveseq,5000),
log='x', type='l', xlab=expression(Wavelength / m),
ylab=expression('Spectral Radiance' / W*sr^{-1}*m^{-2}*m^{-1}), col='blue')
lines(waveseq, cosplanckLawRadWave(waveseq,4000), col='green')
lines(waveseq, cosplanckLawRadWave(waveseq,3000), col='red')
legend('topright', legend=c('3000K','4000K','5000K'), col=c('red','green','blue'), lty=1)

#CMB now:

plot(10^seq(9,12,by=0.01), cosplanckLawRadFreq(10^seq(9,12,by=0.01)),
log='x', type='l', xlab=expression(Frequency / Hz),
ylab=expression('Spectral Radiance' / W*sr^{-1}*m^{-2}*Hz^{-1}))
abline(v=cosplanckPeakFreq(),lty=2)

plot(10^seq(-4,-1,by=0.01), cosplanckLawRadWave(10^seq(-4,-1,by=0.01)),
log='x', type='l', xlab=expression(Wavelength / m),
ylab=expression('Spectral Radiance' / W*sr^{-1}*m^{-2}*m^{-1}))
abline(v=cosplanckPeakWave(),lty=2)

#CMB at surface of last scattering:

TempLastScat=cosplanckCMBTemp(1100) #Note this is still much cooler than our Sun!

plot(10^seq(12,15,by=0.01), cosplanckLawRadFreq(10^seq(12,15,by=0.01),TempLastScat),
log='x', type='l', xlab=expression(Frequency / Hz),
ylab=expression('Spectral Radiance' / W*sr^{-1}*m^{-2}*Hz^{-1}))
abline(v=cosplanckPeakFreq(TempLastScat),lty=2)

plot(10^seq(-7,-4,by=0.01), cosplanckLawRadWave(10^seq(-7,-4,by=0.01),TempLastScat),
log='x', type='l', xlab=expression(Wavelength / m),
ylab=expression('Spectral Radiance' / W*sr^{-1}*m^{-2}*m^{-1}))
abline(v=cosplanckPeakWave(TempLastScat),lty=2)

#Exact number of photons produced by black body:

cosplanckLawRadPhotN()

#We can get pretty much the correct answer through direct integration, i.e.:

integrate(cosplanckLawRadFreqN,1e8,1e12)
integrate(cosplanckLawRadWaveN,1e-4,1e-1)

#Stefan-Boltzmann Law:

cosplanckSBLawRad_sr()

#We can get (almost, some rounding is off) the same answer by multiplying
#the total number of photons produced by a black body per metre squared per
#second per steradian by the average photon energy:

cosplanckLawRadPhotEnAv()*cosplanckLawRadPhotN()

Exact angular area calculator

Description

This function takes a survey geometry defined by RA (long) and Dec (latitude) limits and calculates the exact angular area covered.

Usage

skyarea(long = c(129, 141), lat = c(-2, 3), inunit = "deg", outunit = "deg2", sep=":")

Arguments

Value

Two value vector. First value is the sky area covered in units of outunit (named area), second value is the fraction of the celestial sphere covered by the specified geometry (named areafrac).

Author(s)

Aaron Robotham

See Also

cosvol, hms2deg, dms2deg

Examples

#The GAMA survey areas:
G02area=skyarea(c(30.2,38.8),c(-10.25,-3.72))
G09area=skyarea(c(129,141),c(-2,3))
G12area=skyarea(c(174,186),c(-3,2))
G15area=skyarea(c(211.5,223.5),c(-2,3))
G23area=skyarea(c(338.1,351.9),c(-35,-30))

#Total GAMA survey area:
TotalGAMAarea=G02area+G09area+G12area+G15area+G23area
paste('The GAMA survey area is',round(TotalGAMAarea['area'],2),'sq. deg.')

#Future TACs note: this is less than 1% of the sky ;-)
paste('The GAMA survey area is',round(TotalGAMAarea['areafrac']*100,2),'% of the sky')

Tan Gnomonic and Sine Orthographic Projection System WCS Solver Functions

Description

Converts RA/Dec (degrees) to x/y (pixels) position using the Tan Gnomonic or Sine Orthographic projection systems, and vice-versa. Translations adapted from: http://mathworld.wolfram.com/GnomonicProjection.html and http://mathworld.wolfram.com/OrthographicProjection.html.

Usage

radec2xy(RA, Dec, header, CRVAL1 = 0, CRVAL2 = 0, CRPIX1 = 0, CRPIX2 = 0, CD1_1 = 1,
CD1_2 = 0, CD2_1 = 0, CD2_2 = 1, CTYPE1 = 'RA--TAN', CTYPE2 = 'DEC--TAN')
xy2radec(x, y, header, CRVAL1 = 0, CRVAL2 = 0, CRPIX1 = 0, CRPIX2 = 0, CD1_1 = 1,
CD1_2 = 0, CD2_1 = 0, CD2_2 = 1, CTYPE1 = 'RA--TAN', CTYPE2 = 'DEC--TAN')

Arguments

Details

These functions encode the standard FITS Tan Gnomonic and Sine Orthographic projection systems for solving an image WCS (covering most moden imaging and radio data). They do not deal with higher order polynomial distortion terms.

Value

Author(s)

Aaron Robotham

References

http://mathworld.wolfram.com/GnomonicProjection.html http://mathworld.wolfram.com/OrthographicProjection.html

See Also

deg2dms, deg2hms, dms2deg, hms2deg

Examples

#A simple example:

radec2xy(10, 20)
xy2radec(radec2xy(10, 20))
xy2radec(radec2xy(10, 20, CTYPE1='RA--SIN', CTYPE2='DEC--SIN'),
CTYPE1='RA--SIN',CTYPE2='DEC--SIN')

#A more complicated example, where we transform and rotate large amounts:

exdata_start=expand.grid(1:10,21:30)
plot(exdata_start)
exradec=radec2xy(exdata_start, CRVAL1=20, CRPIX1=100, CRVAL2=30, CRPIX2=130, CD1_1=0.1,
CD1_2=-0.05, CD2_1=0.05, CD2_2=0.1)
plot(exradec)
exdata_end=xy2radec(exradec, CRVAL1=20, CRPIX1=100, CRVAL2=30, CRPIX2=130, CD1_1=0.1,
CD1_2=-0.05, CD2_1=0.05, CD2_2=0.1)
plot(exdata_start,cex=2)
points(exdata_end,col='red')

#The residuals should be very small (in the noice of double precision arithmetic):

plot(density(exdata_start[,1]-exdata_end[,1]))
lines(density(exdata_start[,2]-exdata_end[,2]),col='red')

Transforms 3D spherical coordinates to cartesian coordinates

Description

Transforms 3D spherical coordinates to cartesian coordinates. The user can choose to input the spherical coordinates in degrees or radians.

Usage

sph2car(long, lat, radius = 1, deg = TRUE)

Arguments

Details

This is a low level function that is used for plot transformations.

Value

A data.frame is returned containing the columns x, y and z.

Author(s)

Aaron Robotham

See Also

coordmatch

Examples

print(sph2car(45,0,sqrt(2),deg=TRUE))