Home > reduc > pixdat.m

pixdat

PURPOSE ^

--------------------------------------------------------------------

SYNOPSIS ^

function d = pixdat(d,nside);

DESCRIPTION ^

--------------------------------------------------------------------
d=pixdat(d,nside) 
Script to pixellate to healpix a set of data read in for C-BASS
d is the C-BASS data structure and nside the healpix resolution
Requires the data structure d, the healpix resolution parameter nside and
adds a data column d.antenna0.receiver.pix to the structure with the pixel
number of each point
Charles Copley

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function d = pixdat(d,nside);
0002 %--------------------------------------------------------------------
0003 %d=pixdat(d,nside)
0004 %Script to pixellate to healpix a set of data read in for C-BASS
0005 %d is the C-BASS data structure and nside the healpix resolution
0006 %Requires the data structure d, the healpix resolution parameter nside and
0007 %adds a data column d.antenna0.receiver.pix to the structure with the pixel
0008 %number of each point
0009 %Charles Copley
0010 
0011 %Read in the java wrapper for the healpix library (check it out at h =
0012 %gov.fnal.eag.healpix.PixTools http://home.fnal.gov/~kuropat/HEALPIX/doc/gov/fnal/eag/healpix/PixTools.html#ang2pix_nest%28long,%20double,%20double%29)
0013 
0014 
0015 %Healpix Pixellator
0016 %
0017 % input  : the C-BASS data archive data stored in the d structure
0018             %nside the size of the healpix Nside parameter
0019 % output : Healpix Pixel Numbers
0020 % Tested : matlab 7.8.0
0021 %     By : Charles Copley       May 2010
0022 %--------------------------------------------------------------------
0023 
0024 disp('pixdat:: Calculating the Healpix Pixels in the Data')
0025 h = gov.fnal.eag.healpix.PixTools;
0026 %print('Number of pixels\n')
0027 h.Nside2Npix(nside);
0028 equa = d.antenna0.servo.equa;
0029 long=equa(:,1);
0030 %correct declination into 0->pi (0 at the North pole) ->see
0031 %http://healpix.jpl.nasa.gov/html/java/healpix/core/Healpix.html
0032 colat=(pi()./2-(equa(:,2)));
0033 
0034 pix=zeros(1,length(long));
0035 for i=1:length(long)
0036     pix(i) = h.ang2pix_nest(nside,colat(i),long(i));
0037 end
0038 i=0;
0039 d.antenna0.receiver.pix=pix';
0040 
0041 
0042 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
0043 function OutCoo=horiz_coo(InCoo,JD,TopoPos,Direction);
0044 %--------------------------------------------------------------------
0045 % horiz_coo function        Horizontal coordinates conversion
0046 %                       Converting from equatorial coordinates
0047 %                       to Horizontal coordinates and visa
0048 %                       versa.
0049 % input  : - Two columns matrix of coordinates,
0050 %            (Long & Lat) | (Az & Alt) in radians
0051 %          - vector of JDs + UT fraction of day,
0052 %            if scalar value is given, then it duplicated
0053 %            for all coordinates.
0054 %          - Geodetic Coordinates, east long & north lat in radians
0055 %            if scalar value is given then it is duplicate for
0056 %            all coordinates.
0057 %          - Direction,
0058 %            'h' - from equatorial to horizontal (default).
0059 %            'e' - from horizontal to equatorial.
0060 % output : - two column matrix of output coordinates.
0061 % Tested : matlab 5.3
0062 %     By : Eran O. Ofek           August 1999
0063 %    URL : http://wise-obs.tau.ac.il/~eran/matlab.html
0064 %--------------------------------------------------------------------
0065 if (nargin==3),
0066    Direction = 'h';
0067 elseif (nargin==4),
0068    % no default
0069 else
0070    error('Illigal number of input arguments');
0071 end
0072 
0073 N = length(InCoo(:,1));
0074 if (length(JD)==1),
0075    JD = JD.*ones(N,1);
0076 end
0077 
0078 if (length(TopoPos(:,1))==1),
0079    TopoPos = ones(N,1)*TopoPos;
0080 end
0081 
0082 
0083 % Don't convert Geodetic latitude to Geocentric.
0084 %GeodLat = TopoPos(:,2);
0085 %GeocLatTemp = geod2geoc([TopoPos(:,1),GeodLat,zeros(N,1)],'WGS84');
0086 %GeocLat     = GeocLatTemp(:,2);
0087 %TopoPos(:,2) = GeocLat;
0088 
0089 
0090 % calculating Local Mean Sidereal Time
0091 LST=lst(JD,TopoPos(:,1),'m');
0092 
0093 
0094 if (Direction=='h'),
0095    % convert equatorial to horizontal
0096 
0097    % calculate the Hour Angle
0098    HA  = LST.*2.*pi - InCoo(:,1);
0099    Dec = InCoo(:,2);
0100    Lat = TopoPos(:,2);
0101 
0102    SinAlt = sin(Dec).*sin(Lat) + cos(Dec).*cos(HA).*cos(Lat);
0103    CosAlt = sqrt(1-SinAlt.*SinAlt);
0104 
0105    SinAz  = (-cos(Dec).*sin(HA))./CosAlt;
0106    CosAz  = (sin(Dec).*cos(Lat) - cos(Dec).*cos(HA).*sin(Lat))./CosAlt;
0107 
0108    Az     = atan2(SinAz, CosAz);
0109    Alt    = asin(SinAlt);
0110 
0111    I      = find(Az<0);
0112    Az(I)  = 2.*pi+Az(I);
0113 
0114    OutCoo = [Az, Alt];
0115 elseif (Direction=='e'),
0116    Az     = InCoo(:,1);
0117    Alt    = InCoo(:,2);
0118    Lat = TopoPos(:,2);
0119 
0120    SinDec = sin(Alt).*sin(Lat) + cos(Alt).*cos(Az).*cos(Lat);
0121    CosDec = sqrt(1 - SinDec.*SinDec);
0122 
0123    SinHA  = (-cos(Alt).*sin(Az))./CosDec;
0124    CosHA  = (sin(Alt).*cos(Lat) - cos(Alt).*cos(Az).*sin(Lat))./CosDec;
0125    HA     = atan2(SinHA, CosHA);
0126 
0127    RA     = LST.*2.*pi - HA;
0128    Dec    = asin(SinDec);
0129 
0130    % converting to range [0,1)
0131    RA     = 2.*pi.*(RA./(2.*pi) - floor(RA./(2.*pi)));
0132 
0133    OutCoo = [RA, Dec];
0134 else
0135    error('pixdat:: Illegal type of conversion');
0136 end

Generated on Sun 14-Jun-2015 17:12:45 by m2html © 2005