-------------------------------------------------------------------- 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
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