def calc_rate(redshift):
"""Calculate the rate, kpc/px."""
# Init
# Hubble constant at z = 0
H0 = 71.0
# Omega0, total matter density
Om0 = 0.27
# Cosmo
cosmo = FlatLambdaCDM(H0=H0, Om0=Om0)
# Angular diameter distance, [Mpc]
dA = cosmo.angular_diameter_distance(redshift)
# Resolution of Chandra
rate_px = 0.492 * au.arcsec # [arcmin/px]
rate_px_rad = rate_px.to(au.rad)
# rate
rate = dA.value * rate_px_rad.value # [Mpc/px]
return rate
python类rad()的实例源码
test_visibility_operations.py 文件源码
项目：algorithm-reference-library
作者: SKA-ScienceDataProcessor
项目源码
文件源码
阅读 25
收藏 0
点赞 0
评论 0
def test_phase_rotation(self):
self.vis = create_visibility(self.lowcore, self.times, self.frequency,
channel_bandwidth=self.channel_bandwidth,
phasecentre=self.phasecentre, weight=1.0,
polarisation_frame=PolarisationFrame("stokesIQUV"))
self.vismodel = predict_skycomponent_visibility(self.vis, self.comp)
# Predict visibilities with new phase centre independently
ha_diff = -(self.compabsdirection.ra - self.phasecentre.ra).to(u.rad).value
vispred = create_visibility(self.lowcore, self.times + ha_diff, self.frequency,
channel_bandwidth=self.channel_bandwidth,
phasecentre=self.compabsdirection, weight=1.0,
polarisation_frame=PolarisationFrame("stokesIQUV"))
vismodel2 = predict_skycomponent_visibility(vispred, self.comp)
# Should yield the same results as rotation
rotatedvis = phaserotate_visibility(self.vismodel, newphasecentre=self.compabsdirection, tangent=False)
assert_allclose(rotatedvis.vis, vismodel2.vis, rtol=1e-7)
assert_allclose(rotatedvis.uvw, vismodel2.uvw, rtol=1e-7)
def NES_healpixels(nside=set_default_nside(), width=15, dec_min=0., fill_gap=True):
"""
Define the North Ecliptic Spur region. Return a healpix map with NES pixels as 1.
"""
ra, dec = ra_dec_hp_map(nside=nside)
result = np.zeros(ra.size)
coord = SkyCoord(ra=ra*u.rad, dec=dec*u.rad)
eclip_lat = coord.barycentrictrueecliptic.lat.radian
good = np.where((np.abs(eclip_lat) <= np.radians(width)) & (dec > dec_min))
result[good] += 1
if fill_gap:
good = np.where((dec > np.radians(dec_min)) & (ra < np.radians(180)) &
(dec < np.radians(width)))
result[good] = 1
return result
def galactic_plane_healpixels(nside=set_default_nside(), center_width=10., end_width=4.,
gal_long1=70., gal_long2=290.):
"""
Define the Galactic Plane region. Return a healpix map with GP pixels as 1.
"""
ra, dec = ra_dec_hp_map(nside=nside)
result = np.zeros(ra.size)
coord = SkyCoord(ra=ra*u.rad, dec=dec*u.rad)
g_long, g_lat = coord.galactic.l.radian, coord.galactic.b.radian
good = np.where((g_long < np.radians(gal_long1)) & (np.abs(g_lat) < np.radians(center_width)))
result[good] += 1
good = np.where((g_long > np.radians(gal_long2)) & (np.abs(g_lat) < np.radians(center_width)))
result[good] += 1
# Add tapers
slope = -(np.radians(center_width)-np.radians(end_width))/(np.radians(gal_long1))
lat_limit = slope*g_long+np.radians(center_width)
outside = np.where((g_long < np.radians(gal_long1)) & (np.abs(g_lat) > np.abs(lat_limit)))
result[outside] = 0
slope = (np.radians(center_width)-np.radians(end_width))/(np.radians(360. - gal_long2))
b = np.radians(center_width)-np.radians(360.)*slope
lat_limit = slope*g_long+b
outside = np.where((g_long > np.radians(gal_long2)) & (np.abs(g_lat) > np.abs(lat_limit)))
result[outside] = 0
return result
def is_NES(nside=64, width=15, dec_min=0., fill_gap=True):
"""
Define the North Ecliptic Spur region. Return a healpix map with NES pixels as true.
"""
ra, dec = ra_dec_hp_map(nside=nside)
result = np.zeros(ra.size)
coord = SkyCoord(ra=ra*u.rad, dec=dec*u.rad)
eclip_lat = coord.barycentrictrueecliptic.lat.radian
good = np.where((np.abs(eclip_lat) <= np.radians(width)) & (dec > dec_min))
result[good] += 1
if fill_gap:
good = np.where((dec > np.radians(dec_min)) & (ra < np.radians(180)) &
(dec < np.radians(width)))
result[good] = 1
NES_indx = (result==1)
return NES_indx
def rho_as_angle(rho, eph):
"""Projected linear distance to angular distance.
Parameters
----------
rho : `~astropy.units.Quantity`
Projected distance in units of length.
eph : dictionary-like or `~sbpy.data.Ephem`
Ephemerides; requires geocentric distance as `delta`.
Returns
-------
rho_l : `~astropy.units.Quantity`
"""
if rho.unit.is_equivalent(u.m):
rho_a = np.arctan(rho / eph['delta'].to(u.m))
else:
assert rho.unit.is_equivalent(u.rad)
rho_a = rho
return rho_a
def rho_as_length(rho, eph):
"""Angular distance to projected linear distance.
Parameters
----------
rho : `~astropy.units.Quantity`
Projected distance in units of angle.
eph : dictionary-like or `~sbpy.data.Ephem`
Ephemerides; requires geocentric distance as `delta`.
Returns
-------
rho_l : `~astropy.units.Quantity`
"""
if rho.unit.is_equivalent(u.rad):
rho_l = eph['delta'].to(u.m) * np.tan(rho)
else:
assert rho.unit.is_equivalent(u.m)
rho_l = rho
return rho_l
testing_support.py 文件源码
项目：algorithm-reference-library
作者: SKA-ScienceDataProcessor
项目源码
文件源码
阅读 29
收藏 0
点赞 0
评论 0
def create_configuration_from_file(antfile: str, name: str = None, location: EarthLocation = None,
mount: str = 'altaz',
names: str = "%d", frame: str = 'local',
diameter=35.0,
meta: dict = None,
rmax=None, **kwargs) -> Configuration:
""" Define from a file
:param names:
:param antfile: Antenna file name
:param name: Name of array e.g. 'LOWBD2'
:param location:
:param mount: mount type: 'altaz', 'xy'
:param frame: 'local' | 'global'
:param diameter: Effective diameter of station or antenna
:param meta: Any meta info
:return: Configuration
"""
antxyz = numpy.genfromtxt(antfile, delimiter=",")
assert antxyz.shape[1] == 3, ("Antenna array has wrong shape %s" % antxyz.shape)
if frame == 'local':
latitude = location.geodetic[1].to(u.rad).value
antxyz = xyz_at_latitude(antxyz, latitude)
if rmax is not None:
lantxyz = antxyz - numpy.average(antxyz, axis=0)
r = numpy.sqrt(lantxyz[:, 0] ** 2 + lantxyz[:, 1] ** 2 + lantxyz[:, 2] ** 2)
antxyz = antxyz[r < rmax]
log.debug('create_configuration_from_file: Maximum radius %.1f m includes %d antennas/stations' %
(rmax, antxyz.shape[0]))
nants = antxyz.shape[0]
anames = [names % ant for ant in range(nants)]
mounts = numpy.repeat(mount, nants)
fc = Configuration(location=location, names=anames, mount=mounts, xyz=antxyz, frame=frame,
diameter=diameter)
return fc
def kpc_per_arcsec(self, z):
"""
Calculate the transversal length (unit: kpc) corresponding to
1 arcsec at the *angular diameter distance* of z.
"""
dist_kpc = self.angular_diameter_distance(z, unit="kpc")
return dist_kpc * au.arcsec.to(au.rad)
def kpc_per_pix(self, z):
"""
Calculate the transversal length (unit: kpc) corresponding to
1 ACIS pixel (i.e., 0.492 arcsec) at the *angular diameter distance*
of z.
"""
pix = ACIS.pixel2arcsec * au.arcsec
dist_kpc = self.angular_diameter_distance(z, unit="kpc")
return dist_kpc * pix.to(au.rad).value
def find_contours(data, segments, sigma_level):
"""
This function ...
:param data:
:param segments:
:param sigma_level:
:return:
"""
# Initialize a list for the contours
contours = []
# Get the segment properties
# Since there is only one segment in the source.mask (the center segment), the props
# list contains only one entry (one galaxy)
properties_list = source_properties(data, segments)
for properties in properties_list:
# Obtain the position, orientation and extent
position = Position(properties.xcentroid.value, properties.ycentroid.value)
a = properties.semimajor_axis_sigma.value * sigma_level
b = properties.semiminor_axis_sigma.value * sigma_level
angle = properties.orientation.value # in radians
angle = Angle(angle, u.rad)
radius = Extent(a, b)
meta = {"text": str(properties.label)}
# Create the contour
contours.append(Ellipse(position, radius, angle, meta=meta))
# Return the contours
return contours
# -----------------------------------------------------------------
def find_contour(box, mask, sigma_level):
"""
This function ...
:param box:
:param mask:
:param sigma_level:
:return:
"""
props = source_properties(box, mask)
#tbl = properties_table(props)
x_shift = box.x_min
y_shift = box.y_min
# Since there is only one segment in the self.source.mask (the center segment), the props
# list contains only one entry (one galaxy)
if len(props) == 0: return None
properties = props[0]
# Obtain the position, orientation and extent
position = Position(properties.xcentroid.value + x_shift, properties.ycentroid.value + y_shift)
a = properties.semimajor_axis_sigma.value * sigma_level
b = properties.semiminor_axis_sigma.value * sigma_level
angle = properties.orientation.value # in radians
angle = Angle(angle, u.rad)
radius = Extent(a, b)
# Create and return the elliptical contour
return Ellipse(position, radius, angle)
# -----------------------------------------------------------------
def find_contours(data, segments, sigma_level):
"""
This function ...
:param data:
:param segments:
:param sigma_level:
:return:
"""
# Initialize a list for the contours
contours = []
# Get the segment properties
# Since there is only one segment in the source.mask (the center segment), the props
# list contains only one entry (one galaxy)
properties_list = source_properties(data, segments)
for properties in properties_list:
# Obtain the position, orientation and extent
position = Position(properties.xcentroid.value, properties.ycentroid.value)
a = properties.semimajor_axis_sigma.value * sigma_level
b = properties.semiminor_axis_sigma.value * sigma_level
angle = properties.orientation.value # in radians
angle = Angle(angle, u.rad)
radius = Extent(a, b)
meta = {"text": str(properties.label)}
# Create the contour
contours.append(Ellipse(position, radius, angle, meta=meta))
# Return the contours
return contours
# -----------------------------------------------------------------
def find_contour(box, mask, sigma_level):
"""
This function ...
:param box:
:param mask:
:param sigma_level:
:return:
"""
props = source_properties(box, mask)
#tbl = properties_table(props)
x_shift = box.x_min
y_shift = box.y_min
# Since there is only one segment in the self.source.mask (the center segment), the props
# list contains only one entry (one galaxy)
if len(props) == 0: return None
properties = props[0]
# Obtain the position, orientation and extent
position = Position(properties.xcentroid.value + x_shift, properties.ycentroid.value + y_shift)
a = properties.semimajor_axis_sigma.value * sigma_level
b = properties.semiminor_axis_sigma.value * sigma_level
angle = properties.orientation.value # in radians
angle = Angle(angle, u.rad)
radius = Extent(a, b)
# Create and return the elliptical contour
return Ellipse(position, radius, angle)
# -----------------------------------------------------------------
def calculate(self):
ephem_location = ephem.Observer()
ephem_location.lat = self.location.latitude.to(u.rad) / u.rad
ephem_location.lon = self.location.longitude.to(u.rad) / u.rad
ephem_location.elevation = self.location.height / u.meter
ephem_location.date = ephem.Date(self.time.datetime)
if self.data is None:
self.alt = Latitude([], unit=u.deg)
self.az = Longitude([], unit=u.deg)
self.names = Column([], dtype=np.str)
self.vmag = Column([])
else:
ra = Longitude((self.data['RAh'], self.data['RAm'], self.data['RAs']), u.h)
dec = Latitude((np.core.defchararray.add(self.data['DE-'], self.data['DEd'].astype(str)).astype(int), self.data['DEm'], self.data['DEs']), u.deg)
c = SkyCoord(ra, dec, frame='icrs')
altaz = c.transform_to(AltAz(obstime=self.time, location=self.location))
self.alt = altaz.alt
self.az = altaz.az
self.names = self.data['Name']
self.vmag = self.data['Vmag']
for ephemeris in self.ephemerides:
ephemeris.compute(ephem_location)
self.vmag = self.vmag.insert(0, ephemeris.mag)
self.alt = self.alt.insert(0, (ephemeris.alt.znorm * u.rad).to(u.deg))
self.az = self.az.insert(0, (ephemeris.az * u.rad).to(u.deg))
self.names = self.names.insert(0, ephemeris.name)
return self.names, self.vmag, self.alt, self.az
def selectSrc(self):
"""
Filter a given source, running gtselect
"""
# Do we have to deal with a FITS file or an ASCII list of FITS file ?
allskyext = os.path.splitext(self.allsky)[1]
if allskyext in [".fit", ".fits"]:
fermi.filter['infile'] = self.allsky
else:
fermi.filter['infile'] = '@%s' % self.allsky
if self.daily:
outfile = self.workDir + '/' + str(self.src) + '_daily.fits'
else:
outfile = self.workDir + '/' + str(self.src) + '.fits'
fermi.filter['outfile'] = outfile
# If outfile already exists, we don't do anything
if os.path.isfile(outfile):
return True
fermi.filter['ra'] = self.ra
fermi.filter['dec'] = self.dec
fermi.filter['rad'] = self.roi
fermi.filter['emin'] = self.emin
fermi.filter['emax'] = self.emax
fermi.filter['tmin'] = self.tstart
fermi.filter['tmax'] = self.tstop
fermi.filter['zmax'] = self.zmax
fermi.filter['evclass'] = 128
logging.info('Running gtselect')
fermi.filter.run()
def exposure(self, gamma=None):
"""
Compute exposure on source src, to add a flux column for the photometric light curve.
Warning: the input file is modified in place, with an additional exposure column added to the file !
"""
if self.daily:
infile = self.workDir + '/' + str(self.src) + '_daily_lc.fits'
srcmdl = self.workDir + '/' + str(self.src) + '_daily.xml'
else:
infile = self.workDir + '/' + str(self.src) + '_lc.fits'
srcmdl = self.workDir + '/' + str(self.src) + '.xml'
# If infile already contains an EXPOSURE column, we don't do anything
hdu = fits.open(infile)
if hdu[1].header.get('TTYPE5') == 'EXPOSURE':
return True
scfile = self.spacecraft
irfs = 'P8R2_SOURCE_V6'
rad = str(self.roi)
options = 'infile=' + infile + ' scfile=' + scfile + ' irfs=' + irfs + ' rad=' + rad
if self.fglName is not None:
target = self.fglName.replace('3FGLJ', '3FGL J')
logging.debug('exposure: target=%s', target)
options += ' srcmdl=' + srcmdl + ' target="' + target + '"'
else:
options += ' srcmdl="none" specin=' + str(gamma)
cmd = 'time -p ' + self.fermiDir + '/bin/gtexposure ' + options
logging.info('Running gtexposure')
os.system(cmd)
def stupidFast_altAz2RaDec(alt, az, lat, lon, mjd):
"""
Convert alt, az to RA, Dec without taking into account abberation, precesion, diffraction, ect.
Parameters
----------
alt : numpy.array
Altitude, same length as `ra` and `dec`. Radians.
az : numpy.array
Azimuth, same length as `ra` and `dec`. Must be same length as `alt`. Radians.
lat : float
Latitude of the observatory in radians.
lon : float
Longitude of the observatory in radians.
mjd : float
Modified Julian Date.
Returns
-------
ra : array_like
RA, in radians.
dec : array_like
Dec, in radians.
"""
lmst, last = calcLmstLast(mjd, lon)
lmst = lmst/12.*np.pi # convert to rad
sindec = np.sin(lat)*np.sin(alt) + np.cos(lat)*np.cos(alt)*np.cos(az)
sindec = inrange(sindec)
dec = np.arcsin(sindec)
ha = np.arctan2(-np.sin(az)*np.cos(alt), -np.cos(az)*np.sin(lat)*np.cos(alt)+np.sin(alt)*np.cos(lat))
ra = (lmst-ha)
raneg = np.where(ra < 0)
ra[raneg] = ra[raneg] + 2.*np.pi
raover = np.where(ra > 2.*np.pi)
ra[raover] -= 2.*np.pi
return ra, dec
def is_GP(nside=64, center_width=10., end_width=4.,
gal_long1=70., gal_long2=290.):
"""
Define the Galactic Plane region. Return a healpix map with GP pixels as true.
"""
ra, dec = ra_dec_hp_map(nside=nside)
result = np.zeros(ra.size)
coord = SkyCoord(ra=ra*u.rad, dec=dec*u.rad)
g_long, g_lat = coord.galactic.l.radian, coord.galactic.b.radian
good = np.where((g_long < np.radians(gal_long1)) & (np.abs(g_lat) < np.radians(center_width)))
result[good] += 1
good = np.where((g_long > np.radians(gal_long2)) & (np.abs(g_lat) < np.radians(center_width)))
result[good] += 1
# Add tapers
slope = -(np.radians(center_width)-np.radians(end_width))/(np.radians(gal_long1))
lat_limit = slope*g_long+np.radians(center_width)
outside = np.where((g_long < np.radians(gal_long1)) & (np.abs(g_lat) > np.abs(lat_limit)))
result[outside] = 0
slope = (np.radians(center_width)-np.radians(end_width))/(np.radians(360. - gal_long2))
b = np.radians(center_width)-np.radians(360.)*slope
lat_limit = slope*g_long+b
outside = np.where((g_long > np.radians(gal_long2)) & (np.abs(g_lat) > np.abs(lat_limit)))
result[outside] = 0
GP_inx =(result==1)
return GP_inx
def xarray(shape, yx=[0, 0], th=None):
"""2D array of x values.
Parameters
----------
shape : tuple of int
The shape of the resulting array, (y, x).
yx : tuple of float
Offset the array to align with this y, x center.
th : Quanitity, optional
Place the x-axis along this position angle, measured
counterclockwise from the original x-axis.
Returns
-------
x : ndarray
An array of x values.
Examples
--------
>>> from sbpy.imageanalysis.utils import xarray
>>> x = xarray((10, 10))
>>> x[0, 3]
3
"""
import numpy as np
import astropy.units as u
y, x = np.indices(shape)[-2:]
y = y - yx[0]
x = x - yx[1]
if th is not None:
x = x * np.cos(th.to(u.rad).value) + y * np.sin(th.to(u.rad).value)
return x
def yarray(shape, yx=[0, 0], th=None):
"""2D array of y values.
Parameters
----------
shape : tuple of int
The shape of the resulting array, (y, x).
yx : tuple of float
Offset the array to align with this y, x center.
th : Quantity, optional
Place the y-axis along this position angle, measured
counterclockwise from the original y-axis.
Returns
-------
y : ndarray
An array of y values.
>>> from sbpy.imageanalysis.utils import yarray
>>> y = yarray((10, 10))
>>> y[3, 0]
3
"""
import numpy as np
import astropy.units as u
y, x = np.indices(shape)[-2:]
y = y - yx[0]
x = x - yx[1]
if th is not None:
y = -x * np.sin(th.to(u.rad).value) + y * np.cos(th.to(u.rad).value)
return y
base.py 文件源码
项目：algorithm-reference-library
作者: SKA-ScienceDataProcessor
项目源码
文件源码
阅读 24
收藏 0
点赞 0
评论 0
def phaserotate_visibility(vis: Visibility, newphasecentre: SkyCoord, tangent=True, inverse=False) -> Visibility:
"""
Phase rotate from the current phase centre to a new phase centre
If tangent is False the uvw are recomputed and the visibility phasecentre is updated.
Otherwise only the visibility phases are adjusted
:param vis: Visibility to be rotated
:param newphasecentre:
:param tangent: Stay on the same tangent plane? (True)
:param inverse: Actually do the opposite
:return: Visibility
"""
assert isinstance(vis, Visibility), "vis is not a Visibility: %r" % vis
l, m, n = skycoord_to_lmn(newphasecentre, vis.phasecentre)
# No significant change?
if numpy.abs(n) > 1e-15:
# Make a new copy
newvis = copy_visibility(vis)
phasor = simulate_point(newvis.uvw, l, m)
nvis, npol = vis.vis.shape
# TODO: Speed up (broadcast rules not obvious to me)
if inverse:
for i in range(nvis):
for pol in range(npol):
newvis.data['vis'][i, pol] *= phasor[i]
else:
for i in range(nvis):
for pol in range(npol):
newvis.data['vis'][i, pol] *= numpy.conj(phasor[i])
# To rotate UVW, rotate into the global XYZ coordinate system and back. We have the option of
# staying on the tangent plane or not. If we stay on the tangent then the raster will
# join smoothly at the edges. If we change the tangent then we will have to reproject to get
# the results on the same image, in which case overlaps or gaps are difficult to deal with.
if not tangent:
if inverse:
xyz = uvw_to_xyz(vis.data['uvw'], ha=-newvis.phasecentre.ra.rad, dec=newvis.phasecentre.dec.rad)
newvis.data['uvw'][...] = \
xyz_to_uvw(xyz, ha=-newphasecentre.ra.rad, dec=newphasecentre.dec.rad)[...]
else:
# This is the original (non-inverse) code
xyz = uvw_to_xyz(newvis.data['uvw'], ha=-newvis.phasecentre.ra.rad, dec=newvis.phasecentre.dec.rad)
newvis.data['uvw'][...] = xyz_to_uvw(xyz, ha=-newphasecentre.ra.rad, dec=newphasecentre.dec.rad)[
...]
newvis.phasecentre = newphasecentre
return newvis
else:
return vis
def mergeGTIfiles(self):
"""
Merge multiple GTI files when mergelongterm is True.
Use gtselect.
Assume the current workDir is longTerm/merged.
"""
# Create list of GTI files
if not self.daily:
listname = self.workDir + '/' + self.src + '_gti.list'
else:
listname = self.workDir + '/' + self.src + '_daily_gti.list'
filelist = open(listname, 'w')
list = []
if not self.daily:
for file in glob.glob(self.workDir + '/../20????/' + self.src + '_gti.fits'):
list.append(file)
else:
for file in glob.glob(self.workDir + '/../20????/' + self.src + '_daily_gti.fits'):
list.append(file)
# Sort the list of GTI files
list = sorted(list)
for item in list:
filelist.write(item + '\n')
filelist.close()
fermi.filter['infile'] = '@' + listname
if not self.daily:
outfile = self.workDir + '/' + str(self.src) + '_gti.fits'
else:
outfile = self.workDir + '/' + str(self.src) + '_daily_gti.fits'
fermi.filter['outfile'] = outfile
# If outfile already exists, we re-create it
if os.path.isfile(outfile):
os.remove(outfile)
fermi.filter['ra'] = self.ra
fermi.filter['dec'] = self.dec
fermi.filter['rad'] = self.roi
fermi.filter['emin'] = self.emin
fermi.filter['emax'] = self.emax
fermi.filter['tmin'] = self.tstart
fermi.filter['tmax'] = self.tstop
fermi.filter['zmax'] = self.zmax
fermi.filter['evclass'] = 128
logging.info('Running gtselect')
fermi.filter.run()
def empty_observation():
"""
Return a numpy array that could be a handy observation record
XXX: Should this really be "empty visit"? Should we have "visits" made
up of multple "observations" to support multi-exposure time visits?
XXX-Could add a bool flag for "observed". Then easy to track all proposed
observations. Could also add an mjd_min, mjd_max for when an observation should be observed.
That way we could drop things into the queue for DD fields.
XXX--might be nice to add a generic "sched_note" str field, to record any metadata that
would be useful to the scheduler once it's observed. and/or observationID.
Returns
-------
numpy array
Notes
-----
The numpy fields have the following structure
RA : float
The Right Acension of the observation (center of the field) (Radians)
dec : float
Declination of the observation (Radians)
mjd : float
Modified Julian Date at the start of the observation (time shutter opens)
exptime : float
Total exposure time of the visit (seconds)
filter : str
The filter used. Should be one of u, g, r, i, z, y.
rotSkyPos : float
The rotation angle of the camera relative to the sky E of N (Radians)
nexp : int
Number of exposures in the visit.
airmass : float
Airmass at the center of the field
FWHMeff : float
The effective seeing FWHM at the center of the field. (arcsec)
skybrightness : float
The surface brightness of the sky background at the center of the
field. (mag/sq arcsec)
night : int
The night number of the observation (days)
"""
names = ['RA', 'dec', 'mjd', 'exptime', 'filter', 'rotSkyPos', 'nexp',
'airmass', 'FWHMeff', 'FWHM_geometric', 'skybrightness', 'night', 'slewtime', 'fivesigmadepth',
'alt', 'az', 'clouds', 'moonAlt', 'sunAlt', 'note']
# units of rad, rad, days, seconds, string, radians (E of N?)
types = [float, float, float, float, '|U1', float, int, float, float, float, float, int, float, float,
float, float, float, float, float, '|U40']
result = np.zeros(1, dtype=list(zip(names, types)))
return result
def stupidFast_RaDec2AltAz(ra, dec, lat, lon, mjd, lmst=None):
"""
Convert Ra,Dec to Altitude and Azimuth.
Coordinate transformation is killing performance. Just use simple equations to speed it up
and ignore abberation, precesion, nutation, nutrition, etc.
Parameters
----------
ra : array_like
RA, in radians.
dec : array_like
Dec, in radians. Must be same length as `ra`.
lat : float
Latitude of the observatory in radians.
lon : float
Longitude of the observatory in radians.
mjd : float
Modified Julian Date.
Returns
-------
alt : numpy.array
Altitude, same length as `ra` and `dec`. Radians.
az : numpy.array
Azimuth, same length as `ra` and `dec`. Radians.
"""
if lmst is None:
lmst, last = calcLmstLast(mjd, lon)
lmst = lmst/12.*np.pi # convert to rad
ha = lmst-ra
sindec = np.sin(dec)
sinlat = np.sin(lat)
coslat = np.cos(lat)
sinalt = sindec*sinlat+np.cos(dec)*coslat*np.cos(ha)
sinalt = inrange(sinalt)
alt = np.arcsin(sinalt)
cosaz = (sindec-np.sin(alt)*sinlat)/(np.cos(alt)*coslat)
cosaz = inrange(cosaz)
az = np.arccos(cosaz)
signflip = np.where(np.sin(ha) > 0)
az[signflip] = 2.*np.pi-az[signflip]
return alt, az
