pixdatC

--------------------------------------------------------------------
d=pixdatC(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

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