s2let_hpx_sampling_ring

 s2let_hpx_sampling_ring 
 Compute the Healpix sampling nodes
 Default usage :

   [thetas, phis] = s2let_hpx_sampling_ring(nside)

 nside is the input Healpix resolution,
 thetas contains the colatitudes of the nodes,
 phis contains the longitudes of the nodes.

 S2LET package to perform Wavelets transform on the Sphere.
 Copyright (C) 2012-2015  Boris Leistedt & Jason McEwen
 See LICENSE.txt for license details

0001 function [thetas, phis] = s2let_hpx_sampling_ring(nside)
0002 
0003 % s2let_hpx_sampling_ring
0004 % Compute the Healpix sampling nodes
0005 % Default usage :
0006 %
0007 %   [thetas, phis] = s2let_hpx_sampling_ring(nside)
0008 %
0009 % nside is the input Healpix resolution,
0010 % thetas contains the colatitudes of the nodes,
0011 % phis contains the longitudes of the nodes.
0012 %
0013 % S2LET package to perform Wavelets transform on the Sphere.
0014 % Copyright (C) 2012-2015  Boris Leistedt & Jason McEwen
0015 % See LICENSE.txt for license details
0016 
0017 %     npix = 12 * nside^2;
0018 %     thetas = zeros(npix,1);
0019 %     phis = zeros(npix,1);
0020 %     for i = 1:npix
0021 %        [t,p] = healpix_pix2ang_ring(i, nside);
0022 %        thetas(i) = t;
0023 %        phis(i) = p;
0024 %     end
0025 
0026     npix = 12 * nside^2;
0027     ipix = 1:npix;
0028     nl2 = 2 * nside;
0029     ncap = 2  * nside * (nside - 1);  % points in each polar cap, =0 for nside =1
0030 
0031     thetas = zeros(npix,1);
0032     phis = zeros(npix,1);
0033     
0034     ind = (ipix <= ncap); % North polar cap
0035         hip   = ipix(ind) .* 0.5;
0036         fihip = fix(hip);
0037         iring = fix( sqrt( hip - sqrt(fihip) ) ) + 1 ;% counted from North pole
0038         iphi  = ipix(ind) - 2*iring.*(iring - 1) ;
0039 
0040         thetas(ind) = acos( 1.0 - iring.^2.0 ./ (3.0*nside^2.0) );
0041         phis(ind)   = (iphi - 0.5) .* pi ./ (2.0.*iring);
0042 
0043     ind = (ipix <= nl2*(5*nside+1) & ipix > ncap); % Equatorial region
0044 
0045        ip    = ipix(ind) - ncap - 1;
0046        nl4   = 4*nside;
0047        iring = floor( ip / nl4 ) + nside; % counted from North pole
0048        iphi  = mod(ip,nl4) + 1;
0049 
0050        fodd  = 0.5 * (1 + mod(iring+nside,2)) ; % 1 if iring+nside is odd, 1/2 otherwise
0051        thetas(ind) = acos( (nl2 - iring) ./ (1.5*nside) );
0052        phis(ind)   = (iphi - fodd) * pi /(2.0*nside);
0053 
0054     ind = (ipix > nl2*(5*nside+1) & ipix > ncap); % South polar cap
0055 
0056        ip    = npix - ipix(ind) + 1;
0057        hip   = ip*0.5;
0058        fihip = round(hip);
0059        iring = floor( sqrt( hip - sqrt(fihip) ) ) + 1 ; % counted from South pole
0060        iphi  = 4*iring + 1 - (ip - 2*iring.*(iring-1));
0061 
0062        thetas(ind) = acos( -1.0 + iring.^2 ./ (3.0*nside^2) );
0063        phis(ind)   = (iphi - 0.5) * pi ./ (2.0*iring);
0064 
0065 end
0066 
0067 function [theta, phi] = healpix_pix2ang_ring(ipix, nside)
0068 
0069     npix = 12 * nside^2;
0070     nl2 = 2 * nside;
0071     ncap = 2  * nside * (nside - 1);  % points in each polar cap, =0 for nside =1
0072 
0073     if (ipix <= ncap)
0074         hip   = ipix * 0.5;
0075         fihip = fix(hip);
0076         iring = fix( sqrt( hip - sqrt(fihip) ) ) + 1 ;% counted from North pole
0077         iphi  = ipix - 2*iring*(iring - 1) ;
0078 
0079         theta = acos( 1.0 - iring^2.0 / (3.0*nside^2.0) );
0080         phi   = (iphi - 0.5) * pi/(2.0*iring);
0081 
0082     elseif (ipix <= nl2*(5*nside+1)) % Equatorial region ------
0083 
0084        ip    = ipix - ncap - 1;
0085        nl4   = 4*nside;
0086        iring = floor( ip / nl4 ) + nside; % counted from North pole
0087        iphi  = mod(ip,nl4) + 1;
0088 
0089        fodd  = 0.5 * (1 + mod(iring+nside,2)) ; % 1 if iring+nside is odd, 1/2 otherwise
0090        theta = acos( (nl2 - iring) / (1.5*nside) );
0091        phi   = (iphi - fodd) * pi /(2.0*nside);
0092 
0093     else % South Polar cap -----------------------------------
0094 
0095        ip    = npix - ipix + 1;
0096        hip   = ip*0.5;
0097        fihip = round(hip);
0098        iring = floor( sqrt( hip - sqrt(fihip) ) ) + 1 ; % counted from South pole
0099        iphi  = 4*iring + 1 - (ip - 2*iring*(iring-1));
0100 
0101        theta = acos( -1.0 + iring^2 / (3.0*nside^2) );
0102        phi   = (iphi - 0.5) * pi/(2.0*iring);
0103 
0104     end
0105 
0106 end