Index: /issm/trunk-jpl/externalpackages/mealpix/help/alm2pix_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/alm2pix_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/alm2pix_help.m (revision 25367) @@ -0,0 +1,33 @@ +%% alm2pix +% Evaluate function, expressed as spherical harmonic expansion +% coefficients, on HEALPix pixel centers + +%% Syntax +% v = alm2pix( nSide, alm, 'Param1', Value1, 'Param2', Value2, ...) + +%% Input Arguments +% nSide HEALPix resolution parameter (power of 2) +% alm Spherical harmonic expansion coefficients +% +% Param Value +% 'lmax' max order of harmonic to calculate (default floor2*nSide/3) +% 'nest' (true | {false}) + +%% Return Arguments +% v Function values at pixels in either ring or nest indexing + +%% Example +% Shows the gibbs-like error at the poles +nSide = 2^4; +alm = pix2alm(ones(1,12*nSide^2)); +v = alm2pix(nSide,alm); +hp3d(v-1); + +%% See also +% pix2alm + +%% Requires +% ylm + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/alm2spec_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/alm2spec_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/alm2spec_help.m (revision 25367) @@ -0,0 +1,39 @@ +%% alm2spec +% Evaluate angular power spectral density + +%% Syntax +% cL = alm2spec(alm) +% cL = alm2spec(alm, lMax) + +%% Input Arguments +% alm real valued coefficients +% lMax max order of harmonic to calculate (default 20) + +%% Return Arguments +% cL angular power spectrum + +%% Description +% The angular power spectrum of the function f, +% +% $$f(\theta,\phi) = \sum a_{\ell m}Y_{\ell m}(\theta,\phi)$$ +% +% are the coefficients +% +% $$c_{\ell} = \frac{\sum_{m=-\ell}^\ell |a_{\ell m}|^2}{2*\ell+1}$$ + +%% Example +% Find power spectrum coefficients of dummy data set +nSide = 2^4; +nPix = nSide2nPix(nSide); +pix = 3*(ylm(nSide,4,2) + ylm(nSide,4,1)) + sqrt(17)*ylm(nSide,8,-4) ... + + rand(1,nPix); +lMax = 10; +alm = pix2alm(pix,lMax); +cL = alm2spec(alm,lMax); +stem(0:lMax,cL); + +%% See also +% pix2alm, alm2pix + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/ang2pix_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/ang2pix_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/ang2pix_help.m (revision 25367) @@ -0,0 +1,40 @@ +%% ang2pix +% Convert spherical coordinate location to HEALPix pixel number + +%% Syntax +% pix = ang2pix(nSide, v); +% pix = ang2pix(nSide, v, nest); + +%% Input Arguments +% nSide HEALPix resolution parameter +% v (theta, phi) pairs in radians as either [2,n] or as cellarray +% of [2,1] +% nest (optional) nested pixel number scheme (true | {false}) + +%% Return Arguments +% pix array of pix numbers: +% v on input output +% iscell(v) size(v) cell array of pixel numbers +% isnumeric(v) size(v,1) numeric array of pixel numbers +% +% Note that Matlab numbers pixels from 1 + +%% Example +vTheta = acos(2*rand(3,5)-1); +vPhi = 2*pi*rand(3,5); +vc = arrayfun(@(x,y)([x;y]),vTheta,vPhi,'UniformOutput',false); +va = reshape(cell2mat(vc),[2,3,5]); +nSide = 8; +size(vc) +size(va) +pix = ang2pix(nSide,vc) +pix = ang2pix(nSide,va) + +%% See also +% vec2pix, pix2ang, pix2vec + +%% Requires +% ring2nest + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/ang2vec_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/ang2vec_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/ang2vec_help.m (revision 25367) @@ -0,0 +1,26 @@ +%% ang2vec +% Convert from spherical to cartesian coordinates + +%% Syntax +% xyz = ang2vec(tp) + +%% Input Arguments +% tp cell of [2,1] angular locations (theta; phi) + +%% Return Arguments +% xyz size(tp) celll array of [3,1] vectors + +%% Example +tp = [acos(2*rand(1,12)-1); 2*pi*rand(1,12)]; +tp = mat2cell(tp,2,ones(12,1)); +tp = reshape(tp,3,4); +xyz = ang2vec(tp); +xyz{2,3}' +tp23 = tp{2,3}; +[sin(tp23(1))*cos(tp23(2)),sin(tp23(1))*sin(tp23(2)),cos(tp23(1))] + +%% See also +% vec2ang + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/angDist_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/angDist_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/angDist_help.m (revision 25367) @@ -0,0 +1,66 @@ +%% angDist +% Compute angular distance between pixel centers or points on the sphere + +%% Syntax +% rad = angDist(v0,v1,'Param1',Value1,'Param2',Value2,...) + +%% Input Arguments +% v0, v1 (numerical or cell) array of points on sphere +% +% Param Value +% 'nSide' HEALPix resolution parameter (integer power 2) +% 'nest' nested indexing (true | {false}) + +%% Return Arguments +% rad array of angular distances (radians) + +%% Description +% angDist computes the angular distances between points on the sphere. +% Points may be specified as MEALPix pixel numbers, spherical coordinate +% location vectors, or cartesian coordinate vectors. Cartesian coordinate +% input vectors need not be normalized. +% +% Either or both of v0, v1 may be a scalar; otherwise they must both be +% the same size. When one is a scalar than the distance is computed between +% that point and all the points specified in the complementary argument. +% When pix0 and pix1 are arrays then the distances are returned between the +% corresponding points in the two arrays. +% +% If any location is specified as a pixel number than nSide must be +% specified. The indexing scheme for pixels center locations defaults to +% ring; to specify nested indexing set the 'nest' parameter to true. + +%% v0 or v1 specified as numeric arrays +% If all points in v0 are pixel numbers then v0 may be specified as a +% numeric array (similarly v1, mutatis mundis). then nSide must be specified. +% Pixels numbers may be specified in ring indexing or nested indexing. For +% nested indexing the nest flag must be specified as true. + +%% pix0 or pix1 specified as cell arrays +% Either or both pix0 or pix1 may be specified as a cell array, in which +% case each cell specifies a location on the sphere as either a pixel +% number, a two-component ([theta; phi]) spherical coordinate vector, or a +% three component ([x;y;z]) cartesian vector. nSide must be specified if +% any location is identified as a pixel center. + +%% Example + % angular distance from pole to equator: + v0 = {[0;0;1]}; % north pole + v1 = {[1;0;0]}; % point on equator + rad = angDist(v0, v1); + + % angular distance from pole to pixels 1, 17, 33 for nSide = 4: + v0 = {[0;0;1]}; + rad = angDist([1, 17, 33], v0, 'nSide', 4) + + % angular distance from pixels 1 to 4, 2 to 5, and 3 to 6 for nSide = 8: + rad = angDist([1, 2, 3], [4, 5, 6],'nSide',8) + +%% See also +% queryDisc + +%% Requires +% pix2vec, ang2vec + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/corners_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/corners_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/corners_help.m (revision 25367) @@ -0,0 +1,37 @@ +%% corners +% Find pixel corners + +%% Syntax +% xyz = corners(nSide,nPix,'Param1',Value1,...); + +%% Input Arguments +% nSide HEALPix resolution parameter +% nPix (optional) pixel list. Default all pixels. +% +% Param Value +% 'nest' nested indexing (true | {false}) + +%% Return Arguments +% xyz size(nPix) cell array of [3,1] cartesian vector pixel corner +% locations + +%% Example: +nSide = 2^2; +% Corners for all 48 nSide=4 pixels in ring-indexed order +xyzR = corners(nSide); +size(xyzR) +xyzR{1} +xyzR{end-1} + +nPix = randi(nSide2nPix(nSide), 3, 4); +% Corners for an assortment of pixels in nested index order +xyzN = corners(nSide,nPix,'nest',true); +% Convert to pixel list to ring indexing +r = xyzR{nest2ring(nSide,nPix(2,3))} +% Corners for same pixels requested via ring indexing +n = xyzN{2,3} +% Compare +all(r == n) + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/helpfuncalpha.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/helpfuncalpha.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/helpfuncalpha.m (revision 25367) @@ -0,0 +1,30 @@ +%% Function Alphabetical List +% * - Evaluate spherical harmonic function expansion on HEALPix pixels +% * - valuate angular power spectral density +% * - Convert spherical coordinate location to HEALPix pixel number +% * - Convert from spherical to cartesian coordinates +% * - Computes angular distance on the unit sphere +% * - Find pixel corners +% * - Plots a HEALPix dataset on a 3D sphere +% * - Compute descriptive statistics of a pixel value list +% * - return pixels in ring or ring section +% * - Convert number of pixels to number of facet side subdivisions +% * - Convert resolution parameter to number of pixels +% * - Calculates angular resolution of a HEALPix map in arcsec +% * - Convert MEALPix pixel numbers from nest to ring indexing +% * - Find spherical harmonic decomposition of function on sphere +% * - Find spherical coordinate location of HEALPix pixels +% * - convert HEALPix pixel numbers to (x,y,z) unit vector +% * - Find pixel cartesian face coordinates +% * - Find all pixels whose centers are within a specified disc +% * - Finds minimum base pixel sub-division for required resolution +% * - Convert MEALPix pixel numbers from ring to nest indexing +% * - Find HEALPix ring number corresponding to a given z coordinate +% * - Convert cartesian direction vector to position angle on sphere +% * - Convert cartesian direction vectors to MEALPix pixel numbers +% * - Find pixel number from face and coordinates in face +% * - Return spherical harmonic basis function $Y_L^M$ on the HEALPix sphere +% * - Find HEALPix ring numbers from direction vector z coordinate + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/helpfuncbycat.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/helpfuncbycat.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/helpfuncbycat.m (revision 25367) @@ -0,0 +1,49 @@ +%% Functions by Category +% MEALPix Toolbox +% Version 3.0 +% 21 February 2011 + +%% Coordinate Conversions +% * - Convert from spherical to cartesian coordinates +% * - Convert cartesian direction vector to position angle on sphere + +%% Coordinate-Pixel Conversions +% * - Convert spherical coordinate location to HEALPix pixel number +% * - Find spherical coordinate location of HEALPix pixels +% * - Convert HEALPix pixel numbers to (x,y,z) unit vector +% * - Find HEALPix ring number corresponding to a given z coordinate +% * - Convert cartesian direction vectors to MEALPix pixel numbers +% * - Find HEALPix ring numbers from direction vector z coordinate + +%% Functions on the sphere +% * - Evaluate spherical harmonic function expansion on HEALPix pixels +% * - Computes angular distance on the unit sphere +% * - Find spherical harmonic decomposition of function on sphere +% * - returns spherical harmonic basis $Y_L^M$ on the HEALPix sphere + +%% Map Statistics +% * - valuate angular power spectral density +% * - Compute descriptive statistics of a pixel value list + +%% Pixel Indexing Conversions +% * - Convert MEALPix pixel numbers from nest to ring indexing +% * - Convert MEALPix pixel numbers from ring to nest indexing + +%% Pixel Topology +% * - Find pixel corners +% * - return pixels in ring or ring section +% * - Find pixel cartesian face coordinates +% * - Find all pixels whose centers are within a specified disc +% * - gives the pixel number ipix (NESTED) of ix, iy and face_num + +%% Plotting +% * - Plots a HEALPix dataset on a 3D sphere + +%% Resolution Relations +% * - Convert number of pixels to number of facet side subdivisions +% * - Convert resolution parameter to number of pixels +% * - Calculates angular resolution of a HEALPix map in arcsec +% * - Finds minimum base pixel sub-division for required resolution + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/helptoc.xml =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/helptoc.xml (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/helptoc.xml (revision 25367) @@ -0,0 +1,119 @@ + + + + + + MEALPix Toolbox + + + + + + + + + + + + Getting Started with the MEALPix Toolbox + Overview + System Requirements + Installation + MEALPix for the HEALPix User + Reporting problems + References + Notices + + + + + + Function Reference + + Coordinate Conversions + + ang2vec + vec2ang + + Coordinate-Pixel Conversions + ang2pix + pix2ang + pix2vec + ring2z + vec2pix + z2ring + + Functions on the sphere + alm2pix + angDist + pix2alm + ylm + + Map Statistics + alm2spec + hpStat + + Pixel Indexing Conversions + + nest2ring + ring2nest + + + Pixel Topology + corners + inRing + pix2xy + queryDisc + + xy2pix + + Plotting + + hp3d + + + Resolution relations + nPix2nSide + nSide2nPix + nSide2res + res2nSide + + + + + + + + + + MEALPix Toolbox Web Site + + + Index: /issm/trunk-jpl/externalpackages/mealpix/help/hp3d_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/hp3d_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/hp3d_help.m (revision 25367) @@ -0,0 +1,31 @@ +%% hp3d +% Plot a HEALPix dataset on a 3D sphere + +%% Syntax +% h = hp3d(vPix,'Param1',Value1,'Param2',Value2,...); + +%% Input Arguments +% vPix vector of values at HEALPix pixels +% +% Param Value +% 'nest' vPix in nest index order (true | {false}) + +%% Return Arguments +% h (optional) array of patch handles for the sphere + +%% Description +% Let $f$ be a (real-valued) function on the sphere and vPix be a vector of +% function values at pixel $k$, $k=1\ldots12n_{\textrm{side}}^2$. +% hp3d(vPix,'nest',nest) plots the function on a sphere. The parameter nest +% indicates whether the values vPix are given in nested index order (true) +% or ring index order (false and default). + +%% Example +% plot nested pixel numbers on a sphere +hp3d(1:768,'nest',true); + +%% Requires +% corners, nPix2nSide + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/hpStat_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/hpStat_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/hpStat_help.m (revision 25367) @@ -0,0 +1,30 @@ +%% hpStat +% Compute descriptive statistics of a pixel value list + +%% Syntax +% stats = hpStat(v) + +%% Input Arguments +% v array of HEALPix pixel values + +%% Output Arguments +% stats struct +% field contents +% ntot total number of data points +% nvalid number n of valid data points +% mind minimum of valid data +% maxd maximum of valid data +% average average of valid points +% absdev absolute deviation +% var variance +% stddev standard deviation +% skew skewness factor +% kurt kurtosis factor + +%% Example +% display stats of Y_3^1 +v = ylm(32,3,1,'real',true); +stats = hpStat(abs(v)) + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/alm2pix_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/alm2pix_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/alm2pix_help.html (revision 25367) @@ -0,0 +1,106 @@ + + + + + alm2pix

alm2pix

Evaluate function, expressed as spherical harmonic expansion coefficients, on HEALPix pixel centers

Contents

Syntax

v = alm2pix( nSide, alm, 'Param1', Value1, 'Param2', Value2, ...)

Input Arguments

nSide HEALPix resolution parameter (power of 2) alm Spherical harmonic expansion coefficients

Param Value 'lmax' max order of harmonic to calculate (default floor2*nSide/3) 'nest' (true | {false})

Return Arguments

v       Function values at pixels in either ring or nest indexing

Example

Shows the gibbs-like error at the poles

nSide = 2^4;
+alm = pix2alm(ones(1,12*nSide^2));
+v = alm2pix(nSide,alm);
+hp3d(v-1);
+

See also

pix2alm

Requires

ylm

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/alm2spec_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/alm2spec_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/alm2spec_help.html (revision 25367) @@ -0,0 +1,117 @@ + + + + + alm2spec

alm2spec

Evaluate angular power spectral density

Contents

Syntax

cL = alm2spec(alm)
+cL = alm2spec(alm, lMax)

Input Arguments

alm real valued coefficients lMax max order of harmonic to calculate (default 20)

Return Arguments

cL      angular power spectrum

Description

The angular power spectrum of the function f,

$$f(\theta,\phi) = \sum a_{\ell m}Y_{\ell m}(\theta,\phi)$$

are the coefficients

$$c_{\ell} = \frac{\sum_{m=-\ell}^\ell |a_{\ell m}|^2}{2*\ell+1}$$

Example

Find power spectrum coefficients of dummy data set

nSide = 2^4;
+nPix = nSide2nPix(nSide);
+pix = 3*(ylm(nSide,4,2) + ylm(nSide,4,1)) + sqrt(17)*ylm(nSide,8,-4) ...
+  + rand(1,nPix);
+lMax = 10;
+alm = pix2alm(pix,lMax);
+cL = alm2spec(alm,lMax);
+stem(0:lMax,cL);
+

See also

pix2alm, alm2pix

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/ang2pix_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/ang2pix_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/ang2pix_help.html (revision 25367) @@ -0,0 +1,149 @@ + + + + + ang2pix

ang2pix

Convert spherical coordinate location to HEALPix pixel number

Contents

Syntax

pix = ang2pix(nSide, v);
+pix = ang2pix(nSide, v, nest);

Input Arguments

nSide      HEALPix resolution parameter
+v          (theta, phi) pairs in radians as either [2,n] or as cellarray
+           of [2,1]
+nest       (optional) nested pixel number scheme (true | {false})

Return Arguments

pix             array of pix numbers:
+               v on input     output
+               iscell(v)      size(v) cell array of pixel numbers
+               isnumeric(v)   size(v,1) numeric array of pixel numbers

Note that Matlab numbers pixels from 1

Example

vTheta = acos(2*rand(3,5)-1);
+vPhi = 2*pi*rand(3,5);
+vc = arrayfun(@(x,y)([x;y]),vTheta,vPhi,'UniformOutput',false);
+va = reshape(cell2mat(vc),[2,3,5]);
+nSide = 8;
+size(vc)
+size(va)
+pix = ang2pix(nSide,vc)
+pix = ang2pix(nSide,va)
+
+ans =
+
+     3     5
+
+
+ans =
+
+     2     3     5
+
+
+pix =
+
+   710   602   719   648   504
+   527   167   488    78   610
+   598   437   175   542   152
+
+
+pix =
+
+   710   602   719   648   504
+   527   167   488    78   610
+   598   437   175   542   152
+
+

See also

vec2pix, pix2ang, pix2vec

Requires

ring2nest

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/ang2vec_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/ang2vec_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/ang2vec_help.html (revision 25367) @@ -0,0 +1,112 @@ + + + + + ang2vec

ang2vec

Convert from spherical to cartesian coordinates

Contents

Syntax

xyz = ang2vec(tp)

Input Arguments

tp        cell of [2,1] angular locations (theta; phi)

Return Arguments

xyz       size(tp) celll array of [3,1] vectors

Example

tp = [acos(2*rand(1,12)-1); 2*pi*rand(1,12)];
+tp = mat2cell(tp,2,ones(12,1));
+tp = reshape(tp,3,4);
+xyz = ang2vec(tp);
+xyz{2,3}'
+tp23 = tp{2,3};
+[sin(tp23(1))*cos(tp23(2)),sin(tp23(1))*sin(tp23(2)),cos(tp23(1))]
+
+ans =
+
+   -0.1868    0.1304    0.9737
+
+
+ans =
+
+   -0.1868    0.1304    0.9737
+
+

See also

vec2ang

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/angDist_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/angDist_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/angDist_help.html (revision 25367) @@ -0,0 +1,158 @@ + + + + + angDist

angDist

Compute angular distance between pixel centers or points on the sphere

Contents

Syntax

 rad = angDist(v0,v1,'Param1',Value1,'Param2',Value2,...)

Input Arguments

v0, v1      (numerical or cell) array of points on sphere
Param       Value
+'nSide'     HEALPix resolution parameter (integer power 2)
+'nest'      nested indexing (true | {false})

Return Arguments

 rad     array of angular distances (radians)

Description

angDist computes the angular distances between points on the sphere. Points may be specified as MEALPix pixel numbers, spherical coordinate location vectors, or cartesian coordinate vectors. Cartesian coordinate input vectors need not be normalized.

Either or both of v0, v1 may be a scalar; otherwise they must both be the same size. When one is a scalar than the distance is computed between that point and all the points specified in the complementary argument. When pix0 and pix1 are arrays then the distances are returned between the corresponding points in the two arrays.

If any location is specified as a pixel number than nSide must be specified. The indexing scheme for pixels center locations defaults to ring; to specify nested indexing set the 'nest' parameter to true.

v0 or v1 specified as numeric arrays

If all points in v0 are pixel numbers then v0 may be specified as a numeric array (similarly v1, mutatis mundis). then nSide must be specified. Pixels numbers may be specified in ring indexing or nested indexing. For nested indexing the nest flag must be specified as true.

pix0 or pix1 specified as cell arrays

Either or both pix0 or pix1 may be specified as a cell array, in which case each cell specifies a location on the sphere as either a pixel number, a two-component ([theta; phi]) spherical coordinate vector, or a three component ([x;y;z]) cartesian vector. nSide must be specified if any location is identified as a pixel center.

Example

 % angular distance from pole to equator:
+ v0 = {[0;0;1]}; % north pole
+ v1 = {[1;0;0]}; % point on equator
+ rad = angDist(v0, v1);
+
+ % angular distance from pole to pixels 1, 17, 33 for nSide = 4:
+ v0 = {[0;0;1]};
+ rad = angDist([1, 17, 33], v0, 'nSide', 4)
+
+ % angular distance from pixels 1 to 4, 2 to 5, and 3 to 6 for nSide = 8:
+ rad = angDist([1, 2, 3], [4, 5, 6],'nSide',8)
+
+rad =
+
+    0.2045    0.6224    0.8411
+
+
+rad =
+
+    0.1443    0.2609    0.3013
+
+

See also

queryDisc

Requires

pix2vec, ang2vec

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/corners_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/corners_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/corners_help.html (revision 25367) @@ -0,0 +1,163 @@ + + + + + corners

corners

Find pixel corners

Contents

Syntax

xyz = corners(nSide,nPix,'Param1',Value1,...);

Input Arguments

nSide    HEALPix resolution parameter
+nPix     (optional) pixel list. Default all pixels.
Param    Value
+'nest'   nested indexing (true | {false})

Return Arguments

xyz    size(nPix) cell array of [3,1] cartesian vector pixel corner
+       locations

Example:

nSide = 2^2;
+% Corners for all 48 nSide=4 pixels in ring-indexed order
+xyzR = corners(nSide);
+size(xyzR)
+xyzR{1}
+xyzR{end-1}
+
+nPix = randi(nSide2nPix(nSide), 3, 4);
+% Corners for an assortment of pixels in nested index order
+xyzN = corners(nSide,nPix,'nest',true);
+% Convert to pixel list to ring indexing
+r = xyzR{nest2ring(nSide,nPix(2,3))}
+% Corners for same pixels requested via ring indexing
+n = xyzN{2,3}
+% Compare
+all(r == n)
+
+ans =
+
+     1   192
+
+
+ans =
+
+         0    0.2031    0.2826    0.0000
+         0         0    0.2826    0.2031
+    1.0000    0.9792    0.9167    0.9792
+
+
+ans =
+
+   -0.2826   -0.2031         0   -0.0000
+   -0.2826    0.0000         0   -0.2031
+   -0.9167   -0.9792   -1.0000   -0.9792
+
+
+r =
+
+    0.5478    0.7071    0.5478    0.3827
+    0.8198    0.7071    0.8198    0.9239
+    0.1667         0   -0.1667         0
+
+
+n =
+
+    0.5478    0.7071    0.5478    0.3827
+    0.8198    0.7071    0.8198    0.9239
+    0.1667         0   -0.1667         0
+
+
+ans =
+
+     1     1     1     1
+
+

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/helpfuncalpha.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/helpfuncalpha.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/helpfuncalpha.html (revision 25367) @@ -0,0 +1,100 @@ + + + + + Function Alphabetical List

Function Alphabetical List

  • alm2pix - Evaluate spherical harmonic function expansion on HEALPix pixels
  • alm2spec - valuate angular power spectral density
  • ang2pix - Convert spherical coordinate location to HEALPix pixel number
  • ang2vec - Convert from spherical to cartesian coordinates
  • angDist - Computes angular distance on the unit sphere
  • corners - Find pixel corners
  • hp3d - Plots a HEALPix dataset on a 3D sphere
  • hpStat - Compute descriptive statistics of a pixel value list
  • inRing - return pixels in ring or ring section
  • nPix2nSide - Convert number of pixels to number of facet side subdivisions
  • nSide2nPix - Convert resolution parameter to number of pixels
  • nSide2res - Calculates angular resolution of a HEALPix map in arcsec
  • nest2ring - Convert MEALPix pixel numbers from nest to ring indexing
  • pix2alm - Find spherical harmonic decomposition of function on sphere
  • pix2ang - Find spherical coordinate location of HEALPix pixels
  • pix2vec - convert HEALPix pixel numbers to (x,y,z) unit vector
  • pix2xy - Find pixel cartesian face coordinates
  • queryDisc - Find all pixels whose centers are within a specified disc
  • res2nSide - Finds minimum base pixel sub-division for required resolution
  • ring2nest - Convert MEALPix pixel numbers from ring to nest indexing
  • ring2z - Find HEALPix ring number corresponding to a given z coordinate
  • vec2ang - Convert cartesian direction vector to position angle on sphere
  • vec2pix - Convert cartesian direction vectors to MEALPix pixel numbers
  • xy2pix - Find pixel number from face and coordinates in face
  • ylm - Return spherical harmonic basis function $Y_L^M$ on the HEALPix sphere
  • z2ring - Find HEALPix ring numbers from direction vector z coordinate

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/helpfuncbycat.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/helpfuncbycat.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/helpfuncbycat.html (revision 25367) @@ -0,0 +1,119 @@ + + + + + Functions by Category

Functions by Category

MEALPix Toolbox Version 3.0 21 February 2011

Contents

Coordinate Conversions

  • ang2vec - Convert from spherical to cartesian coordinates
  • vec2ang - Convert cartesian direction vector to position angle on sphere

Coordinate-Pixel Conversions

  • ang2pix - Convert spherical coordinate location to HEALPix pixel number
  • pix2ang - Find spherical coordinate location of HEALPix pixels
  • pix2vec - Convert HEALPix pixel numbers to (x,y,z) unit vector
  • ring2z - Find HEALPix ring number corresponding to a given z coordinate
  • vec2pix - Convert cartesian direction vectors to MEALPix pixel numbers
  • z2ring - Find HEALPix ring numbers from direction vector z coordinate

Functions on the sphere

  • alm2pix - Evaluate spherical harmonic function expansion on HEALPix pixels
  • angDist - Computes angular distance on the unit sphere
  • pix2alm - Find spherical harmonic decomposition of function on sphere
  • ylm - returns spherical harmonic basis $Y_L^M$ on the HEALPix sphere

Map Statistics

  • alm2spec - valuate angular power spectral density
  • hpStat - Compute descriptive statistics of a pixel value list

Pixel Indexing Conversions

  • nest2ring - Convert MEALPix pixel numbers from nest to ring indexing
  • ring2nest - Convert MEALPix pixel numbers from ring to nest indexing

Pixel Topology

  • corners - Find pixel corners
  • inRing - return pixels in ring or ring section
  • pix2xy - Find pixel cartesian face coordinates
  • queryDisc - Find all pixels whose centers are within a specified disc
  • xy2pix - gives the pixel number ipix (NESTED) of ix, iy and face_num

Plotting

  • hp3d - Plots a HEALPix dataset on a 3D sphere

Resolution Relations

  • nPix2nSide - Convert number of pixels to number of facet side subdivisions
  • nSide2nPix - Convert resolution parameter to number of pixels
  • nSide2res - Calculates angular resolution of a HEALPix map in arcsec
  • res2nSide - Finds minimum base pixel sub-division for required resolution

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/hp3d_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/hp3d_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/hp3d_help.html (revision 25367) @@ -0,0 +1,102 @@ + + + + + hp3d

hp3d

Plot a HEALPix dataset on a 3D sphere

Contents

Syntax

h = hp3d(vPix,'Param1',Value1,'Param2',Value2,...);

Input Arguments

vPix   vector of values at HEALPix pixels
Param      Value
+'nest'     vPix in nest index order (true | {false})

Return Arguments

h      (optional) array of patch handles for the sphere

Description

Let $f$ be a (real-valued) function on the sphere and vPix be a vector of function values at pixel $k$, $k=1\ldots12n_{\textrm{side}}^2$. hp3d(vPix,'nest',nest) plots the function on a sphere. The parameter nest indicates whether the values vPix are given in nested index order (true) or ring index order (false and default).

Example

plot nested pixel numbers on a sphere

hp3d(1:768,'nest',true);
+

Requires

corners, nPix2nSide

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/hpStat_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/hpStat_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/hpStat_help.html (revision 25367) @@ -0,0 +1,126 @@ + + + + + hpStat

hpStat

Compute descriptive statistics of a pixel value list

Contents

Syntax

stats = hpStat(v)

Input Arguments

v         array of HEALPix pixel values

Output Arguments

stats     struct
+          field     contents
+          ntot      total number of data points
+          nvalid    number n of valid data points
+          mind      minimum of valid data
+          maxd      maximum of valid data
+          average   average of valid points
+          absdev    absolute deviation
+          var       variance
+          stddev    standard deviation
+          skew      skewness factor
+          kurt      kurtosis factor

Example

display stats of Y_3^1

v = ylm(32,3,1,'real',true);
+stats = hpStat(abs(v))
+
+stats = 
+
+      ntot: 12288
+    nValid: 12288
+      mind: 1.0813e-18
+      maxd: 0.6289
+      mean: 0.2272
+    absdev: 0.1421
+       var: 0.0280
+    stddev: 0.1673
+      skew: 0.5111
+      kurt: 2.2410
+
+

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/inRing_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/inRing_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/inRing_help.html (revision 25367) @@ -0,0 +1,189 @@ + + + + + inRing

inRing

Return list of pixels in a ring or ring section

Contents

Syntax

pixList = inRing(nSide, nRing, phi0, dPhi, 'Param1', Value1, ...);

Input Arguments

nSide       HEALPix resolution parameter
+nRing       Ring number (1 <= nRing <= 4*nSide-1)
+phi0, dPhi  (optional) Ring section longitudes [phi-dPhi,phi+dPhi]
Param       Value
+'nest'      return pixels with nested indexing (true | {false})
+'strict'    return only pixels whose center is in longitude range
+            (true | {false})

Return Arguments

pixList     size(nRing) cell array of pixels in ring defined by nRing,
+            phi0, dPhi

Description

nRing may be a numeric array, in which case phi0, dPhi and nest may each
+be either scalar or a size(nRing) array. If strict is false (default)
+then all pixels intersected by longitude range are included in pixList.

Example

% All pixels in ring 2 of nSide = 4 pixelization (ring indexed)
+pix = inRing(4,2);
+pix{:}
+
+% Same, but nested indexing
+pix = inRing(4,2,'nest',true);
+pix{:}
+nest2ring(4,pix{:})
+
+% Pixels in rings [2,4;8,10;12,14] and the longitude band
+% [7*pi/4,9*pi/4]
+nRing = [2,3;5,12;13,14];
+pix = inRing(4,nRing,2*pi,pi/4);
+pix{1}
+pix{end}
+
+% Same, with strict selection
+pix = inRing(4,nRing,2*pi,pi/4,'strict',true);
+pix{1}
+pix{end}
+
+ans =
+
+     5     6     7     8     9    10    11    12
+
+
+ans =
+
+    65    16    15    32    31    48    47    64
+
+
+ans =
+
+   137     1     5     2     7     3     9     4
+
+
+ans =
+
+    12     5     6
+
+
+ans =
+
+   188   181   182
+
+
+ans =
+
+    12     5
+
+
+ans =
+
+   188   181
+
+

See also

ringNum

Requires

ring2z, pix2vec

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_bugreport.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_bugreport.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_bugreport.html (revision 25367) @@ -0,0 +1,79 @@ + + + + + Reporting bugs

Reporting bugs

MEALPix software distributed on an "as is" basis, without warranties or conditions of any kind, and without promise of support. Nevertheless, the Gravitational Wave Astronomy Group makes extensive use of MEALPix and welcomes any and all bug reports via e-mail to our request-tracking system rt-mealpix@gwastro.org. Bugs reports will be investigated and addressed on a best effort basis.

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_gs_top.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_gs_top.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_gs_top.html (revision 25367) @@ -0,0 +1,80 @@ + + + + + Getting Started with MEALPix

Getting Started with MEALPix

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_healpix.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_healpix.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_healpix.html (revision 25367) @@ -0,0 +1,90 @@ + + + + + MEALPix for the HEALPix User

MEALPix for the HEALPix User

Contents

MEALPix function names are derived from their HEALPix counterparts.

MEALPix functions are typically camel-cased (e.g., angDist, or nSide2nPix). Underscores in HEALPix names are dropped in the conversion to MEALPix names.

MEALPix functions accept array arguments

MEALPix functions exploit Matlab(tm)'s vectorization facilities to allow for the processing of arrays of pixels, direction vectors and position angle vectors via a single function call. For example, the MEALPix implementation of angDist accepts mixed arrays of cartesian direction vectors, angular position vectors, and ring or nest indexed pixels for either argument.

MEALPix indexes pixels from 1

HEALPix pixels are indexed from 0. In keeping with the Matlab standard MEALPix pixels are indexed from 1. MEALPix pixels are indexed from 1: i.e., MEALPix pixel number n corresponds to HEALPix pixel number n-1.

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_install.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_install.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_install.html (revision 25367) @@ -0,0 +1,85 @@ + + + + + Installation

Installation

  • Download MEALPix, which is provided as a gzipped tar file. Un-tar the file and place the MEALPix directory in some appropriate location. Hereafter we will refer to the fully-qualified path to MEALPix as $MEALPIX.
  • Add the MEALPix directory tree to your Matlab™ path using, e.g., the Matlab™ command addpath(genpath('$MEALPIX')). You may wish to add this command to your matlab/startup.m file, in which case MEALPix will be available whenever you start-up Matlab™.

The MEALPix Toolbox includes documentation that will appear in the Matlab Documentation browser when Matlab™ starts-up with MEALPix in its path.

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_notices.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_notices.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_notices.html (revision 25367) @@ -0,0 +1,86 @@ + + + + + Notices

Notices

MEALPix is based on HEALPix. HEALPix is distributed under the terms of the GNU General Public License ("GPL") as published by the Free Software Foundation; either version 2 of the GPL, or (at your option) any later version.

This Matlab™ expression of HEALPix is Copyright 2010-2011 by Lee Samuel Finn. It is licensed under the Apache License, Version 2.0 (the "License"); you may not use MEALPix except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, MEALPix software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_overview.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_overview.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_overview.html (revision 25367) @@ -0,0 +1,101 @@ + + + + + Product Overview

Product Overview

MEALPix is a native Matlab implementation of a substantial subset of the HEALPix subroutine library. It supports the fast and accurate statistical and astrophysical analysis of massive full-sky data sets. HEALPix is supported by NASA and ESO, used extensively by the WMAP and PLANCK missions, and is part of the official FITS World Coordinate System.

MEALPix provides fully vectorized functions for

  • Conversion between position angles, cartesian direction vectors, and HEALPix ring and nest pixel indexing schemes;
  • Spherical harmonic resolution and analysis of functions on the sphere;
  • Nearest pixel identification (pixels within a given disc, pixels in a ring, pixels in a HEALPix face).

MEALPix was created from the Fortran90 implementation HEALPix via a combination of hand coding and machine assisted conversion of the HEALPix F90 source code. Some routines were hand-coded based on the documented functionality of the F90 subroutines on which they were based. Other routines were hand-coded via inspection of the HEALPix F90 implementation source. Still other routines were machine converted from F90 to Matlab™ using f2matlab. In this latter case the resulting Matlab™ source was then hand-modified to vectorize the computations and take advantage of Matlab™ data types and other functions.

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_product_page.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_product_page.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_product_page.html (revision 25367) @@ -0,0 +1,100 @@ + + + + + mealpix_product_page

+Gravitational Wave Astronomy Group Logo +

MEALPix Toolbox

+

Contents

The MEALPix Toolbox is a native Matlab™ implementation of the Hierarchical Equal Area isoLatitude Pixelization of the sphere (HEALPix). It is derived from the HEALPix package, originally developed by ESO and currently supported by JPL and ESO.

Functions

Getting Started

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_references.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_references.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_references.html (revision 25367) @@ -0,0 +1,102 @@ + + + + + References

References

K. M. Gorski et al. 2005. HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere. The Astrophysical Journal 622 759

K. M. Gorski et al. 1999. The HEALPix Primer. arXiv:astro-ph/9905275

K. M. Gorski et al. 2002. HEALPix - a Framework for High Resolution, Fast Analysis on the Sphere. In Astronomical Data Analysis Software and Systems XI, Bohlender, D. A. et al. ed. ASP Conference Series Vol 281, pp. 107-111

Mark R. Calabrettal and Boudewijn F. Roukema. 2007. Mapping on the HEALPix grid. Monthly Notices of the Royal Astronomical Society 381:2 pp. 865 - 872.

W. O'Mullane et al. 2001. Splitting the Sky - HTM and HEALPix. In Mining The Sky. Anthony Banday et al., Ed. Anthony ESO Astrophysics Symposia 2001, 638 - 648.

L. Nicastro and Calderone, G. 2008. Indexing Astronomical Database Tables using HTM and HEALPix. In Astronomical Data Analysis Software and Systems (ADASS) XVII. A. W. Argle et. al, ed. Astronomical Society of the Pacific Conference Series vol. 394 pp. 486-490.

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_reqts.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_reqts.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/mealpix_reqts.html (revision 25367) @@ -0,0 +1,93 @@ + + + + + System Requirements

System Requirements

MEALPix is a native Matlab™ implementation of HEALPix: neither a HEALPix installation nor any third-party java or mex files are required to make use of the basic HEALPix functionality.

Contents

Matlab™ requirements

MEALPix was developed and tested using Matlab™ R2010a. It should work with any later version.

MEALPix does not require any Mathworks™ Toolboxes or add-ons. MEALPix will make use of the capabilities of the Parallel Computing Toolbox, should it be present.

Third-party software requirements

MEALPix also includes functions not found in HEALPix. Some of these additional functions require other toolboxes:

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/nPix2nSide_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/nPix2nSide_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/nPix2nSide_help.html (revision 25367) @@ -0,0 +1,98 @@ + + + + + nPix2nSide

nPix2nSide

Finds number of HEALPix facet side subdivisions from number of pixels on sphere

Contents

Syntax

nSide = nPix2nSide(nPix);

Input Arguments

nPix     number of pixels on sphere

Return Arguments

nSide    sub-divisions of HEALPix facet side

Example

nPix = 12*[1:4].^2;
+nSide = nPix2nSide(nPix)
+
+nSide =
+
+     1     2     3     4
+
+

See also

nSide2nPix

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/nSide2nPix_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/nSide2nPix_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/nSide2nPix_help.html (revision 25367) @@ -0,0 +1,101 @@ + + + + + nSide2nPix

nSide2nPix

Convert HEALPix resolution parameter to number of pixels on sphere

Contents

Syntax

nPix = nSide2nPix(nSide)

Input Arguments

nSide    sub-divisions of facet sides. May be an array.

Return Arguments

nPix     size(nSide) number of pixels on sphere corresponding to nSide

Description

HEALPix divides the sphere into 12 basic facets. Each facet is further divided into nSide sub-divisions for a total of 12*nSide^2 pixels.

Example

nSide = 2.^[1:4];
+nPix = nSide2nPix(nSide)
+
+nPix =
+
+          48         192         768        3072
+
+

See also

nPix2nSide

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/nSide2res_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/nSide2res_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/nSide2res_help.html (revision 25367) @@ -0,0 +1,97 @@ + + + + + nSide2res

nSide2res

Determine angular resolution of a HEALPix map

Contents

Syntax

angres = nSide2res(nSide)

Input Arguments

nSide   HEALPix resolution parameter (may be an array)

Return Arguments

angres  Angular resolution in arcsec

Example

nSide = 2.^(8:12);
+nSide2res(nSide)
+
+ans =
+
+  824.5167  412.2584  206.1292  103.0646   51.5323
+
+

See also

res2nSide

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/nest2ring_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/nest2ring_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/nest2ring_help.html (revision 25367) @@ -0,0 +1,106 @@ + + + + + nest2ring

nest2ring

Convert MEALPix pixel numbers from nest to ring indexing

Contents

Syntax

rPix = nest2ring(nSide,nPix)

Input Arguments

nSide     HEALPix resolution parameters
+nPix      numeric array of MEALPix nest indexed pixel numbers

Return Arguments

rPix      size(nPix) numeric array of ring indexed pixel numbers

Example

nPix = 1:10;
+rPix = nest2ring(2,nPix)
+nPix = ring2nest(2,rPix)
+
+rPix =
+
+    14     6     5     1    16     8     7     2    18    10
+
+
+nPix =
+
+     1     2     3     4     5     6     7     8     9    10
+
+

See also

ring2nest

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/pix2alm_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/pix2alm_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/pix2alm_help.html (revision 25367) @@ -0,0 +1,132 @@ + + + + + pix2alm

pix2alm

Find spherical harmonic decomposition of function on sphere

Contents

Syntax

alm = pix2alm(v)
+alm = pix2alm(v, lmax)

Input Arguments

v       array of pixel values
+lMax    (optional) max order of harmonic to calculate
nPix = numel(v), with nPix = 12*nSide^2 for nSide a power of 2.
+lMax defaults to 2*floor(nSide/3)

Return Arguments

alm     coefficients of spherical harmonic expansion

Description

Let $x_k$ denote the location of pixel $k$ and $v_k$ the function value at $x_k$. Then

$$alm(j) = \frac{4\pi}{N}\sum_{k=0}^{N-1} Y_{LM}^{\dagger}(x_k) v_k$$

where $j=(L+1)^2+M-L$ and $N$ is the number of pixels (12*nSide^2)

Example

estimate alm of dummy data

ns = 2^4;
+np = 12*ns^2;
+v = ylm(ns,1,1) + ylm(ns,2,-2) + ylm(ns,3,0) + rand(1,np)/10;
+lMax = 4;
+alm = pix2alm(v,lMax);
+for L = 0:3
+  fprintf('L = %d: ',L);
+  fprintf('%7.3f ',abs(alm((L+1)^2+(-L:L)-L)));
+  fprintf('\n');
+end
+
L = 0:   0.176 
+L = 1:   0.001   0.002   1.000 
+L = 2:   1.001   0.002   0.001   0.002   0.003 
+L = 3:   0.003   0.003   0.002   0.998   0.002   0.003   0.003 
+

See also

alm2pix

Requires

ylm

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/pix2ang_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/pix2ang_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/pix2ang_help.html (revision 25367) @@ -0,0 +1,135 @@ + + + + + pix2ang

pix2ang

Find spherical coordinate location of HEALPix pixels

Contents

Syntax

tp = pix2ang(nSide,pix,'Param1',Value1,...)

Input Arguments

nSide    HEALPix resolution parameter (positive integer power of 2)
+pix      (optional) array of pixel values to find coordinates of
Param    Values
+'nest'   pixels are given in nested indexing (true | {false})

Return Arguments

tp       cell array of [theta; phi] pixel locations (radians)

Example

Return [theta; phi] postion angles for all pix of nSide=4, nest scheme

nSide = 4;
+tp = pix2ang(nSide,'nest',true);
+size(tp)
+tp{1}
+% Return position angles for a random assortment of pixels
+pix = randi(nSide2nPix(nSide),3,4);
+tp = pix2ang(nSide,pix);
+size(tp)
+tp{1}
+
+ans =
+
+     1   192
+
+
+ans =
+
+    1.4033
+    0.7854
+
+
+ans =
+
+     3     4
+
+
+ans =
+
+    1.5708
+    6.0868
+
+

See also

See also ang2pix, pix2vec, vec2pix

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/pix2vec_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/pix2vec_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/pix2vec_help.html (revision 25367) @@ -0,0 +1,144 @@ + + + + + pix2vec

pix2vec

Convert MEALPix pixel numbers to cartesian location vectors

Contents

Syntax

[x,y,z] = pix2vec(nSide,pix,'Param1',Value1,...)
+xyz = pix2vec(nSide,pix,'Param1',Value1,...)

Input Arguments

nSide       HEALPix resolution parameter
+pix         (optional) numeric array of pixel numbers
Param       Value
+'nest'      indexing scheme is nested (true | {false})

Return Arguments

x,y,z       size(pix) size(pix) numeric array of unit vector coordinate
+            components corresponding to pixels pix
+xyz         size(pix) size(pix) cell array of cartesian unit vectors
+            corresponding to pixels pix

Example

nSide = 2^4;
+nPix = nSide2nPix(nSide);
+pix0 = reshape(1:6,3,2);
+pix1 = reshape(nPix+(-5:0),3,2);
+pix = [pix0,pix1];
+xyz = pix2vec(nSide,pix);
+xyzn = pix2vec(nSide,pix,'nest',true);
+size(xyz)
+[xyz{1},xyz{end}]
+[xyzn{1},xyzn{end}]
+
+ans =
+
+     3     4
+
+
+ans =
+
+    0.0361    0.0361
+    0.0361   -0.0361
+    0.9987   -0.9987
+
+
+ans =
+
+    0.7065    0.7065
+    0.7065   -0.7065
+    0.0417   -0.0417
+
+

See also

ang2pix, pix2ang, vec2pix

Requires

nest2ring

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/pix2xy_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/pix2xy_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/pix2xy_help.html (revision 25367) @@ -0,0 +1,107 @@ + + + + + pix2xy

pix2xy

Find pixel cartesian face coordinates

Contents

Syntax

xy = pix2xy(nSide,nPix)

Input Arguments

nSide     HEALPix resolution parameter
+nPix      nested indexing pixel numbers to find coordinates of

Return Arguments

xy        size(nPix) cell array of pixel coordinates on face

Example

nPix = reshape(1:24,[4,3,2]);
+xy = pix2xy(2,nPix)
+
+xy(:,:,1) = 
+
+    [2x1 double]    [2x1 double]    [2x1 double]
+    [2x1 double]    [2x1 double]    [2x1 double]
+    [2x1 double]    [2x1 double]    [2x1 double]
+    [2x1 double]    [2x1 double]    [2x1 double]
+
+
+xy(:,:,2) = 
+
+    [2x1 double]    [2x1 double]    [2x1 double]
+    [2x1 double]    [2x1 double]    [2x1 double]
+    [2x1 double]    [2x1 double]    [2x1 double]
+    [2x1 double]    [2x1 double]    [2x1 double]
+
+

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/queryDisc_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/queryDisc_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/queryDisc_help.html (revision 25367) @@ -0,0 +1,145 @@ + + + + + queryDisc

queryDisc

Find all pixels whose centers are within a specified disc

Contents

Syntax

pix = queryDisc(nSide,c0,r,'Param1',Value1,...)

Input Arguments

nSide     HEALPix resolution parameter
+c0        disc center. May be specified as a pixel number, spherical
+          coordinate location [theta, phi] in rads, or as a cartesian
+          unit vector [x;y;z].
+r         disc radius (radians)
Param     Value
+'nest'    pixels are returned using nested indexing (true | {false})

Return Arguments

pix      list of pixels in disc

Example

nSide=32;
+% Pixels within pi/10 radians of north pole
+pix = queryDisc(nSide,[0;0;1],pi/10);
+numel(pix)
+nSide2nPix(nSide)*((pi/10)/(4*pi))
+% Pixels within pi/10 radians of south pole
+pix = queryDisc(nSide,[0;0;-1],pi/10);
+numel(pix)
+% Pixels within pi/10 radians of 0 long, lat
+pix = queryDisc(nSide,[1;0;0],pi/10);
+numel(pix)
+
+ans =
+
+   312
+
+
+ans =
+
+  307.2000
+
+
+ans =
+
+   312
+
+
+ans =
+
+   302
+
+

See also

queryPolygon, queryStrip, queryTriangle

Requires

pix2vec, ang2vec, angDist

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/res2nSide_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/res2nSide_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/res2nSide_help.html (revision 25367) @@ -0,0 +1,99 @@ + + + + + res2nSide

res2nSide

Find minimum HEALPix resolution parameter yielding required resolution

Contents

Syntax

nSide = res2nSide(r)

Input Arguments

r        angular resolution (arc seconds)

Return Arguments

nSide    minimum base pixel sub-divisions corresponding to resolution r

Example

r = 100*rand(3,4);
+nSide = res2nSide(r)
+
+nSide =
+
+        8192        4096        4096        4096
+       16384        4096       16384       16384
+        8192        8192        4096       32768
+
+

See also

nSide2res

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/ring2nest_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/ring2nest_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/ring2nest_help.html (revision 25367) @@ -0,0 +1,106 @@ + + + + + ring2nest

ring2nest

Convert MEALPix pixel numbers from ring to nest indexing

Contents

Syntax

nPix = ring2nest(nSide,rPix)

Input Arguments

nSide      HEALPix resolution parameter
+rPix       ring indexed MEALPix pixel numbers

Return Arguments

nPix       next indexed MEALPix pixel numbers

Example

nPix = ring2nest(4,reshape(1:6,2,3))
+rPix = nest2ring(4,nPix)
+
+nPix =
+
+    16    48    15
+    32    64    14
+
+
+rPix =
+
+     1     3     5
+     2     4     6
+
+

See also

nest2ring

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/ring2z_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/ring2z_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/ring2z_help.html (revision 25367) @@ -0,0 +1,108 @@ + + + + + ring2z

ring2z

Find HEALPix ring number corresponding to a given z coordinate

Contents

Syntax

z = ring2z(nSide,nRing)

Input Arguments

nSide    HEALPix resolution parameter
+nRing    numeric array of ring numbers (in range [1,4*nSide-1])

Return Arguments

z        size(nRing) numeric array of pixel center z coordinates

Example

nRing = reshape(1:15,3,5)
+z = ring2z(4,nRing)
+
+nRing =
+
+     1     4     7    10    13
+     2     5     8    11    14
+     3     6     9    12    15
+
+
+z =
+
+    0.9792    0.6667    0.1667   -0.3333   -0.8125
+    0.9167    0.5000         0   -0.5000   -0.9167
+    0.8125    0.3333   -0.1667   -0.6667   -0.9792
+
+

See also

z2ring

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/ringNum_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/ringNum_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/ringNum_help.html (revision 25367) @@ -0,0 +1,145 @@ + + + + + ringNum

ringNum

Find pixel ring numbers

Contents

Syntax

nRing = ringNum(nSide,pix)

Input Arguments

nSide      HEALPix resolution parameter
+pix        (optional) MEALPix pixel number array (ring indexing).
+           Defaults to all pixels

Return Arguments

nRing      size(pix) array of ring numbers

Description

Algorithm Notes

Example

ringNum(2)
+
Undefined function or method 'ringNum' for input arguments of type 'double'.
+
+Error in ==> ringNum_help at 20
+ringNum(2)
+

See also

inRing

Requires

pix2vec, z2ring

References

K. M. Gorski et al. 2005. HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere. The Astrophysical Journal 622 759

K. M. Gorski et al. 1999. The HEALPix Primer. arXiv:astro-ph/9905275

K. M. Gorski et al. 2002. HEALPix - a Framework for High Resolution, Fast Analysis on the Sphere. In Astronomical Data Analysis Software and Systems XI, Bohlender, D. A. et al. ed. ASP Conference Series Vol 281, pp. 107-111

Mark R. Calabrettal and Boudewijn F. Roukema. 2007. Mapping on the HEALPix grid. Monthly Notices of the Royal Astronomical Society 381:2 pp. 865 - 872.

W. O'Mullane et al. 2001. Splitting the Sky - HTM and HEALPix. In Mining The Sky. Anthony Banday et al., Ed. Anthony ESO Astrophysics Symposia 2001, 638 - 648.

L. Nicastro and Calderone, G. 2008. Indexing Astronomical Database Tables using HTM and HEALPix. In Astronomical Data Analysis Software and Systems (ADASS) XVII. A. W. Argle et. al, ed. Astronomical Society of the Pacific Conference Series vol. 394 pp. 486-490.

Notices

Copyright 2010-2011 Lee Samuel Finn Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.nSideWarn(nSide,mfilename);


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/synopsis.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/synopsis.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/synopsis.html (revision 25367) @@ -0,0 +1,127 @@ + + + + + Synopsis

Synopsis

 alm2pix    - Evaluate spherical harmonic function expansion on HEALPix pixels
+ alm2spec   - valuate angular power spectral density
+ ang2pix    - Convert spherical coordinate location to HEALPix pixel number
+ ang2vec    - Convert from spherical to cartesian coordinates
+ angDist    - Computes angular distance on the unit sphere
+ corners    - Find pixel corners
+ hp3d       - Plots a HEALPix dataset on a 3D sphere
+ hpInterp   - Interpolate data on sphere to MEALPix pixels
+ hpStat     - Compute descriptive statistics of a pixel value list
+ inRing     - return pixels in ring or ring section
+ nPix2nSide - Convert number of pixels to number of facet side subdivisions
+ nSide2nPix - Convert resolution parameter to number of pixels
+ nSide2res  - Calculates angular resolution of a HEALPix map in arcsec
+ nest2ring  - Convert MEALPix pixel numbers from nest to ring indexing
+ pix2alm    - Find spherical harmonic decomposition of function on sphere
+ pix2ang    - Find spherical coordinate location of HEALPix pixels
+ pix2vec    - convert HEALPix pixel numbers to (x,y,z) unit vector
+ pix2xy     - PIX2XY
+ queryDisc  - Find all pixels whose centers are within a specified disc
+ res2nSide  - Finds minimum base pixel sub-division for required resolution
+ ring2nest  - Convert MEALPix pixel numbers from ring to nest indexing
+ ring2z     - Find HEALPix ring number corresponding to a given z coordinate
+ vec2ang    - Convert cartesian direction vector to position angle on sphere
+ vec2pix    - Convert cartesian direction vectors to MEALPix pixel numbers
+ xy2pix     - gives the pixel number ipix (NESTED) of ix, iy and face_num
+ ylm        - returns spherical harmonic basis Y_L^M on the HEALPix sphere
+ z2ring     - Find HEALPix ring numbers from direction vector z coordinate

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/template_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/template_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/template_help.html (revision 25367) @@ -0,0 +1,133 @@ + + + + + NAME

NAME

short descr

Contents

Syntax

Input Arguments

Return Arguments

Description

Algorithm Notes

Example

See also

Requires

References

K. M. Gorski et al. 2005. HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere. The Astrophysical Journal 622 759

K. M. Gorski et al. 1999. The HEALPix Primer. arXiv:astro-ph/9905275

K. M. Gorski et al. 2002. HEALPix - a Framework for High Resolution, Fast Analysis on the Sphere. In Astronomical Data Analysis Software and Systems XI, Bohlender, D. A. et al. ed. ASP Conference Series Vol 281, pp. 107-111

Mark R. Calabrettal and Boudewijn F. Roukema. 2007. Mapping on the HEALPix grid. Monthly Notices of the Royal Astronomical Society 381:2 pp. 865 - 872.

W. O'Mullane et al. 2001. Splitting the Sky - HTM and HEALPix. In Mining The Sky. Anthony Banday et al., Ed. Anthony ESO Astrophysics Symposia 2001, 638 - 648.

L. Nicastro and Calderone, G. 2008. Indexing Astronomical Database Tables using HTM and HEALPix. In Astronomical Data Analysis Software and Systems (ADASS) XVII. A. W. Argle et. al, ed. Astronomical Society of the Pacific Conference Series vol. 394 pp. 486-490.

Notices

Copyright 2010-2011 Lee Samuel Finn Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/vec2ang_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/vec2ang_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/vec2ang_help.html (revision 25367) @@ -0,0 +1,135 @@ + + + + + vec2ang

vec2ang

Convert from cartesian direction vector to position angle on sphere

Contents

Syntax

tp = vec2ang(xyz)

Input Arguments

xyz    cell array of [3,1] cartesian coordinate direction vectors

The direction vectors need not be normalized.

Return Arguments

tp     cell array of [2,1] position angles ([theta, phi]) on sphere

Example

Create a [3,4] cell array of position angles and manually convert to cartesian direction vectors. Use vec2ang to convert back to position angles and compare with originals

% angular position vectors
+z = 2*rand(3,4)-1;
+th = acos(z);
+ph = 2*pi*rand(3,4)-pi;
+
+% convert to cartesian coordinate vectors
+x = sqrt(1-z.^2).*cos(ph);
+y = sqrt(1-z.^2).*sin(ph);
+xyz = arrayfun(@(x,y,z)([x,y,z]),x,y,z,'UniformOutput',false);
+
+% use vec2ang to convert cartesian coordinate vectors to angular position
+% vectors
+tp = vec2ang(xyz);
+
+% Compare original vectors with final vectors
+tp0 = arrayfun(@(t,p)([t;p]),th,ph,'UniformOutput',false);
+tf = cellfun(@(x,y)(any(abs(x-y)>4*eps)),tp,tp0)
+
+tf =
+
+     0     0     0     0
+     0     0     0     0
+     0     0     0     0
+
+

See also

ang2vec

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/vec2pix_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/vec2pix_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/vec2pix_help.html (revision 25367) @@ -0,0 +1,133 @@ + + + + + vec2pix

vec2pix

Convert cartesian direction vectors to MEALPix pixel numbers

Contents

Syntax

nPix = vec2pix(nSide,xyz,'Param1',Value1,...);

Input Arguments

nSide     HEALPix resolution parameter
+xyz       cell array of [3,1] cartesian direction vectors
Param     Value
+'nest'    use nested indexing (true | {false})

Return Arguments

nPix      size(xyz) pixel number array

Example

Find pixels corresponding to a [4,5] set of (unnormalized) cartesian vectors:

nSide = 4;
+xyz0 = squeeze(num2cell(2*rand(3,4,5)-1,1));
+nPix = vec2pix(nSide,xyz0)
+
+% Find the normalized vectors corresponding to the pixels
+xyz1 = pix2vec(nSide,nPix);
+
+% Find angular distance from vectors to pixels
+d = angDist(xyz0,xyz1)
+
+nPix =
+
+    28   157     2   112    58
+   120    37    44   140   153
+    82   121    68   127    25
+   165    29    51   140    42
+
+
+d =
+
+    0.1107    0.0752    0.1586    0.0250    0.1124
+    0.0695    0.1483    0.1005    0.1260    0.1347
+    0.0588    0.0572    0.0838    0.1210    0.0761
+    0.1156    0.1455    0.1688    0.0376    0.1256
+
+

See also

ang2pix, pix2ang, pix2vec

Requires

ring2nest

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/xy2pix_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/xy2pix_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/xy2pix_help.html (revision 25367) @@ -0,0 +1,119 @@ + + + + + xy2pix

xy2pix

Find (nest indexed) MEALPix pixel number corresponding to a pixel index on a HEALPix face

Contents

Syntax

ipix = xy2pix(nSide, ix, iy, nFace)

Input Arguments

nside     HEALPix resolution parameter (power of 2)
+ix,iy     integer pixel location on face (1 <= ix, iy <= nSide)
+nFace     HEALPix face  number (1 <= nFace <= 12)

One or both ix and iy may be scalars. If neither is a scalar both must have the same size.

Return Arguments

ipix      NESTED pixel index for given x,y on the face

Description

At lowest resolution HEALPix divides the sphere into 12 faces, each with its own coordinate system. Each face is further subdivided in a regular fashion as the resolution parameter is increased. Given a face number (from 1 to 12), a resolution parameter, and a location on the face XY2PIX identifies the pixel number (in the nested indexing scheme).

Example

nSide = 4;
+nFace = 3;
+ix = [1 2; 1 2];
+iy = [1 1; 2 2];
+nPix = xy2pix(nSide,ix,iy,nFace)
+
+nPix =
+
+    49    50
+    51    52
+
+

See also

pix2xy

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/ylm_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/ylm_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/ylm_help.html (revision 25367) @@ -0,0 +1,103 @@ + + + + + ylm

ylm

Evaluate spherical harmonic basis functions $Y_{LM}$ at pixel centers

Contents

Syntax

v = ylm(nSide,L,M,'Param1',Value1,'Param2',Value2,...);

Input Arguments

nSide    HEALPix resolution parameter (power of 2)
+L        spherical harmonic degree (0 <= L)
+M        spherical harmonic order (-M <= L <= M)
Param    Value
+'real'   return real values (true | {false})
+'nest'   nested indexing flag (true | {false})

Return Arguments

v       spherical harmonics evaluated at pixel centers

Example

Plot Ylm for (L,M) = (4,2) on nSide = 16 HEALPix sphere

pix = ylm(16,4,2);
+hp3d(pix);
+

Requires

pix2ang

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/html/z2ring_help.html =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/html/z2ring_help.html (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/html/z2ring_help.html (revision 25367) @@ -0,0 +1,106 @@ + + + + + z2ring

z2ring

Return HEALPix ring number for a given z

Contents

Syntax

zNum = z2ring(nSide,z)

Input Arguments

nSide      HEALPix resolution parameter
+z          numeric array of z values

Return Arguments

nRing      size(z) numeric array of ring numbers

Description

HEALPix pixels centers are arranged in rings of constant z. Each ring spans a given range of z. z2ring returns the ring numbers that cover the input z values.

Example

nSide = 4;
+zNum = z2ring(nSide,reshape((-7:7)*2/15,3,5))
+
+zNum =
+
+    14    11     9     6     4
+    13    10     8     6     3
+    12    10     7     5     2
+
+

See also

ring2z

Copyright 2010-2011 Lee Samuel Finn. Terms of Use.


+ Published with MATLAB® 7.12

Index: /issm/trunk-jpl/externalpackages/mealpix/help/inRing_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/inRing_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/inRing_help.m (revision 25367) @@ -0,0 +1,57 @@ +%% inRing +% Return list of pixels in a ring or ring section + +%% Syntax +% pixList = inRing(nSide, nRing, phi0, dPhi, 'Param1', Value1, ...); + +%% Input Arguments +% nSide HEALPix resolution parameter +% nRing Ring number (1 <= nRing <= 4*nSide-1) +% phi0, dPhi (optional) Ring section longitudes [phi-dPhi,phi+dPhi] +% +% Param Value +% 'nest' return pixels with nested indexing (true | {false}) +% 'strict' return only pixels whose center is in longitude range +% (true | {false}) + +%% Return Arguments +% pixList size(nRing) cell array of pixels in ring defined by nRing, +% phi0, dPhi + +%% Description +% nRing may be a numeric array, in which case phi0, dPhi and nest may each +% be either scalar or a size(nRing) array. If strict is false (default) +% then all pixels intersected by longitude range are included in pixList. + +%% Example + +% All pixels in ring 2 of nSide = 4 pixelization (ring indexed) +pix = inRing(4,2); +pix{:} + +% Same, but nested indexing +pix = inRing(4,2,'nest',true); +pix{:} +nest2ring(4,pix{:}) + +% Pixels in rings [2,4;8,10;12,14] and the longitude band +% [7*pi/4,9*pi/4] +nRing = [2,3;5,12;13,14]; +pix = inRing(4,nRing,2*pi,pi/4); +pix{1} +pix{end} + +% Same, with strict selection +pix = inRing(4,nRing,2*pi,pi/4,'strict',true); +pix{1} +pix{end} + + +%% See also +% ringNum + +%% Requires +% ring2z, pix2vec + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_bugreport.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_bugreport.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_bugreport.m (revision 25367) @@ -0,0 +1,10 @@ +%% Reporting bugs +% MEALPix software distributed on an "as is" basis, without warranties +% or conditions of any kind, and without promise of support. Nevertheless, +% the makes +% extensive use of MEALPix and welcomes any and all bug reports via e-mail +% to our request-tracking system rt-mealpix@gwastro.org. Reported bugs +% will be investigated and addressed on a best effort basis. + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_gs_top.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_gs_top.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_gs_top.m (revision 25367) @@ -0,0 +1,11 @@ +%% Getting Started with MEALPix +% * +% * +% * +% * +% * +% * +% * + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_healpix.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_healpix.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_healpix.m (revision 25367) @@ -0,0 +1,21 @@ +%% MEALPix for the HEALPix User +%% MEALPix function names are derived from their HEALPix counterparts. +% MEALPix functions are typically camel-cased (e.g., , or ). Underscores in HEALPix +% names are dropped in the conversion to MEALPix names. + +%% MEALPix functions accept array arguments +% MEALPix functions exploit Matlab(tm)'s vectorization facilities to +% allow for the processing of arrays of pixels, direction vectors and +% position angle vectors via a single function call. For example, the +% MEALPix implementation of accepts mixed arrays of +% cartesian direction vectors, angular position vectors, and ring or nest +% indexed pixels for either argument. + +%% MEALPix indexes pixels from 1 +% HEALPix pixels are indexed from 0. In keeping with the Matlab standard +% MEALPix pixels are indexed from 1. MEALPix pixels are indexed from 1: +% i.e., MEALPix pixel number n corresponds to HEALPix pixel number n-1. + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_install.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_install.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_install.m (revision 25367) @@ -0,0 +1,16 @@ +%% Installation +% * Download MEALPix, which is provided as a gzipped tar file. Un-tar the +% file and place the MEALPix directory in some appropriate location. +% Hereafter we will refer to the fully-qualified path to MEALPix as +% |$MEALPIX|. +% +% * Add the MEALPix directory tree to your Matlab(TM) path using, e.g., the +% Matlab(TM) command |addpath(genpath('$MEALPIX'))|. You may wish to add +% this command to your matlab/startup.m file, in which case MEALPix will be +% available whenever you start-up Matlab(TM). +% +% The MEALPix Toolbox includes documentation that will appear in the Matlab +% Documentation browser when Matlab(TM) starts-up with MEALPix in its path. + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_notices.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_notices.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_notices.m (revision 25367) @@ -0,0 +1,17 @@ +%% Notices +% MEALPix is based on HEALPix. HEALPix is distributed under the terms of +% the ("GPL") as published by the Free Software Foundation; +% either version 2 of the GPL, or (at your option) any later version. +% +% This Matlab(TM) expression of HEALPix is Copyright 2010-2011 by Lee +% Samuel Finn. It is licensed under the Apache License, Version 2.0 (the +% "License"); you may not use MEALPix except in compliance with the +% License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. Index: /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_overview.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_overview.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_overview.m (revision 25367) @@ -0,0 +1,32 @@ +%% Product Overview +% MEALPix is a native Matlab implementation of a substantial subset of the +% subroutine library. It supports +% the fast and accurate statistical and astrophysical analysis of massive +% full-sky data sets. HEALPix is supported by NASA and ESO, used +% extensively by the WMAP and PLANCK missions, and is part of the official +% FITS World Coordinate System. +% +% MEALPix provides fully vectorized functions for +% +% * Conversion between position angles, cartesian direction vectors, and +% HEALPix ring and nest pixel indexing schemes; +% +% * Spherical harmonic resolution and analysis of functions on the sphere; +% +% * Nearest pixel identification (pixels within a given disc, pixels in a +% ring, pixels in a HEALPix face). +% +% MEALPix was created from the Fortran90 implementation +% via a combination of hand coding +% and machine assisted conversion of the HEALPix F90 source code. Some +% routines were hand-coded based on the documented functionality of the F90 +% subroutines on which they were based. Other routines were hand-coded via +% inspection of the HEALPix F90 implementation source. Still other routines +% were machine converted from F90 to Matlab(TM) using +% . In +% this latter case the resulting Matlab(TM) source was then hand-modified +% to vectorize the computations and take advantage of Matlab(TM) data types +% and other functions. + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_product_page.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_product_page.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_product_page.m (revision 25367) @@ -0,0 +1,27 @@ +%% +% +% Gravitational Wave Astronomy Group Logo +%

MEALPix Toolbox

+% + +%% +% The MEALPix Toolbox is a native Matlab(TM) implementation of the +% Hierarchical Equal Area isoLatitude Pixelization of the sphere (HEALPix). +% It is derived from the package, +% originally developed by ESO and currently supported by JPL and ESO. + +%% Functions +% * +% * + +%% Getting Started +% * +% * +% * +% * +% * +% * + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_references.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_references.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_references.m (revision 25367) @@ -0,0 +1,32 @@ +%% References +% K. M. Gorski et al. 2005. . The Astrophysical Journal 622 759 +% +% K. M. Gorski et al. 1999. . arXiv:astro-ph/9905275 +% +% K. M. Gorski et al. 2002. +% . +% In Astronomical Data Analysis Software and Systems XI, Bohlender, D. A. +% et al. ed. ASP Conference Series Vol 281, pp. 107-111 +% +% Mark R. Calabrettal and Boudewijn F. Roukema. 2007. +% . Monthly Notices of the Royal Astronomical Society 381:2 +% pp. 865 - 872. +% +% W. O'Mullane et al. 2001. . In Mining The Sky. Anthony Banday +% et al., Ed. Anthony ESO Astrophysics Symposia 2001, 638 - 648. +% +% L. Nicastro and Calderone, G. 2008. +% . In Astronomical Data +% Analysis Software and Systems (ADASS) XVII. A. W. Argle et. al, ed. +% Astronomical Society of the Pacific Conference Series vol. 394 pp. +% 486-490. + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_reqts.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_reqts.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/mealpix_reqts.m (revision 25367) @@ -0,0 +1,23 @@ +%% System Requirements +% MEALPix is a native Matlab(TM) implementation of HEALPix: neither a +% HEALPix installation nor any third-party java or mex files are required +% to make use of the basic HEALPix functionality. + +%% Matlab(TM) requirements +% MEALPix was developed and tested using Matlab(TM) R2010a. It should work +% with any later version. +% +% MEALPix does not require any Mathworks(TM) Toolboxes or add-ons. MEALPix +% will make use of the capabilities of the +% , should it be present. + +%% Third-party software requirements +% MEALPix also includes functions not found in HEALPix. Some of these +% additional functions require other toolboxes: +% +% * requires the +% Toolbox. + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/nPix2nSide_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/nPix2nSide_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/nPix2nSide_help.m (revision 25367) @@ -0,0 +1,22 @@ +%% nPix2nSide +% Finds number of HEALPix facet side subdivisions from number of pixels on +% sphere + +%% Syntax +% nSide = nPix2nSide(nPix); + +%% Input Arguments +% nPix number of pixels on sphere + +%% Return Arguments +% nSide sub-divisions of HEALPix facet side + +%% Example +nPix = 12*[1:4].^2; +nSide = nPix2nSide(nPix) + +%% See also +% nSide2nPix + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/nSide2nPix_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/nSide2nPix_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/nSide2nPix_help.m (revision 25367) @@ -0,0 +1,25 @@ +%% nSide2nPix +% Convert HEALPix resolution parameter to number of pixels on sphere + +%% Syntax +% nPix = nSide2nPix(nSide) + +%% Input Arguments +% nSide sub-divisions of facet sides. May be an array. + +%% Return Arguments +% nPix size(nSide) number of pixels on sphere corresponding to nSide + +%% Description +% HEALPix divides the sphere into 12 basic facets. Each facet is further +% divided into nSide sub-divisions for a total of 12*nSide^2 pixels. + +%% Example +nSide = 2.^[1:4]; +nPix = nSide2nPix(nSide) + +%% See also +% nPix2nSide + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/nSide2res_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/nSide2res_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/nSide2res_help.m (revision 25367) @@ -0,0 +1,21 @@ +%% nSide2res +% Determine angular resolution of a HEALPix map + +%% Syntax +% angres = nSide2res(nSide) + +%% Input Arguments +% nSide HEALPix resolution parameter (may be an array) + +%% Return Arguments +% angres Angular resolution in arcsec + +%% Example +nSide = 2.^(8:12); +nSide2res(nSide) + +%% See also +% res2nSide + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/nest2ring_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/nest2ring_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/nest2ring_help.m (revision 25367) @@ -0,0 +1,23 @@ +%% nest2ring +% Convert MEALPix pixel numbers from nest to ring indexing + +%% Syntax +% rPix = nest2ring(nSide,nPix) + +%% Input Arguments +% nSide HEALPix resolution parameters +% nPix numeric array of MEALPix nest indexed pixel numbers + +%% Return Arguments +% rPix size(nPix) numeric array of ring indexed pixel numbers + +%% Example +nPix = 1:10; +rPix = nest2ring(2,nPix) +nPix = ring2nest(2,rPix) + +%% See also +% ring2nest + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/pix2alm_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/pix2alm_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/pix2alm_help.m (revision 25367) @@ -0,0 +1,46 @@ +%% pix2alm +% Find spherical harmonic decomposition of function on sphere + +%% Syntax +% alm = pix2alm(v) +% alm = pix2alm(v, lmax) + +%% Input Arguments +% v array of pixel values +% lMax (optional) max order of harmonic to calculate +% +% nPix = numel(v), with nPix = 12*nSide^2 for nSide a power of 2. +% lMax defaults to 2*floor(nSide/3) + +%% Return Arguments +% alm coefficients of spherical harmonic expansion + +%% Description +% Let $x_k$ denote the location of pixel $k$ and $v_k$ the function value +% at $x_k$. Then +% +% $$alm(j) = \frac{4\pi}{N}\sum_{k=0}^{N-1} Y_{LM}^{\dagger}(x_k) v_k$$ +% +% where $j=(L+1)^2+M-L$ and $N$ is the number of pixels (12*nSide^2) + +%% Example +% estimate alm of dummy data +ns = 2^4; +np = 12*ns^2; +v = ylm(ns,1,1) + ylm(ns,2,-2) + ylm(ns,3,0) + rand(1,np)/10; +lMax = 4; +alm = pix2alm(v,lMax); +for L = 0:3 + fprintf('L = %d: ',L); + fprintf('%7.3f ',abs(alm((L+1)^2+(-L:L)-L))); + fprintf('\n'); +end + +%% See also +% alm2pix + +%% Requires +% ylm + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/pix2ang_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/pix2ang_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/pix2ang_help.m (revision 25367) @@ -0,0 +1,33 @@ +%% pix2ang +% Find spherical coordinate location of HEALPix pixels + +%% Syntax +% tp = pix2ang(nSide,pix,'Param1',Value1,...) + +%% Input Arguments +% nSide HEALPix resolution parameter (positive integer power of 2) +% pix (optional) array of pixel values to find coordinates of +% +% Param Values +% 'nest' pixels are given in nested indexing (true | {false}) + +%% Return Arguments +% tp cell array of [theta; phi] pixel locations (radians) + +%% Example +% Return [theta; phi] postion angles for all pix of nSide=4, nest scheme +nSide = 4; +tp = pix2ang(nSide,'nest',true); +size(tp) +tp{1} +% Return position angles for a random assortment of pixels +pix = randi(nSide2nPix(nSide),3,4); +tp = pix2ang(nSide,pix); +size(tp) +tp{1} + +%% See also +% See also ang2pix, pix2vec, vec2pix + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/pix2vec_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/pix2vec_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/pix2vec_help.m (revision 25367) @@ -0,0 +1,40 @@ +%% pix2vec +% Convert MEALPix pixel numbers to cartesian location vectors + +%% Syntax +% [x,y,z] = pix2vec(nSide,pix,'Param1',Value1,...) +% xyz = pix2vec(nSide,pix,'Param1',Value1,...) + +%% Input Arguments +% nSide HEALPix resolution parameter +% pix (optional) numeric array of pixel numbers +% +% Param Value +% 'nest' indexing scheme is nested (true | {false}) + +%% Return Arguments +% x,y,z size(pix) size(pix) numeric array of unit vector coordinate +% components corresponding to pixels pix +% xyz size(pix) size(pix) cell array of cartesian unit vectors +% corresponding to pixels pix + +%% Example +nSide = 2^4; +nPix = nSide2nPix(nSide); +pix0 = reshape(1:6,3,2); +pix1 = reshape(nPix+(-5:0),3,2); +pix = [pix0,pix1]; +xyz = pix2vec(nSide,pix); +xyzn = pix2vec(nSide,pix,'nest',true); +size(xyz) +[xyz{1},xyz{end}] +[xyzn{1},xyzn{end}] + +%% See also +% ang2pix, pix2ang, vec2pix + +%% Requires +% nest2ring + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/pix2xy_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/pix2xy_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/pix2xy_help.m (revision 25367) @@ -0,0 +1,19 @@ +%% pix2xy +% Find pixel cartesian face coordinates + +%% Syntax +% xy = pix2xy(nSide,nPix) + +%% Input Arguments +% nSide HEALPix resolution parameter +% nPix nested indexing pixel numbers to find coordinates of + +%% Return Arguments +% xy size(nPix) cell array of pixel coordinates on face + +%% Example +nPix = reshape(1:24,[4,3,2]); +xy = pix2xy(2,nPix) + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/queryDisc_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/queryDisc_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/queryDisc_help.m (revision 25367) @@ -0,0 +1,40 @@ +%% queryDisc +% Find all pixels whose centers are within a specified disc + +%% Syntax +% pix = queryDisc(nSide,c0,r,'Param1',Value1,...) + +%% Input Arguments +% nSide HEALPix resolution parameter +% c0 disc center. May be specified as a pixel number, spherical +% coordinate location [theta, phi] in rads, or as a cartesian +% unit vector [x;y;z]. +% r disc radius (radians) +% +% Param Value +% 'nest' pixels are returned using nested indexing (true | {false}) + +%% Return Arguments +% pix list of pixels in disc + +%% Example +nSide=32; +% Pixels within pi/10 radians of north pole +pix = queryDisc(nSide,[0;0;1],pi/10); +numel(pix) +nSide2nPix(nSide)*((pi/10)/(4*pi)) +% Pixels within pi/10 radians of south pole +pix = queryDisc(nSide,[0;0;-1],pi/10); +numel(pix) +% Pixels within pi/10 radians of 0 long, lat +pix = queryDisc(nSide,[1;0;0],pi/10); +numel(pix) + +%% See also +% queryPolygon, queryStrip, queryTriangle + +%% Requires +% pix2vec, ang2vec, angDist + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/res2nSide_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/res2nSide_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/res2nSide_help.m (revision 25367) @@ -0,0 +1,21 @@ +%% res2nSide +% Find minimum HEALPix resolution parameter yielding required resolution + +%% Syntax +% nSide = res2nSide(r) + +%% Input Arguments +% r angular resolution (arc seconds) + +%% Return Arguments +% nSide minimum base pixel sub-divisions corresponding to resolution r + +%% Example +r = 100*rand(3,4); +nSide = res2nSide(r) + +%% See also +% nSide2res + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/ring2nest_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/ring2nest_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/ring2nest_help.m (revision 25367) @@ -0,0 +1,22 @@ +%% ring2nest +% Convert MEALPix pixel numbers from ring to nest indexing + +%% Syntax +% nPix = ring2nest(nSide,rPix) + +%% Input Arguments +% nSide HEALPix resolution parameter +% rPix ring indexed MEALPix pixel numbers + +%% Return Arguments +% nPix next indexed MEALPix pixel numbers + +%% Example +nPix = ring2nest(4,reshape(1:6,2,3)) +rPix = nest2ring(4,nPix) + +%% See also +% nest2ring + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/ring2z_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/ring2z_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/ring2z_help.m (revision 25367) @@ -0,0 +1,22 @@ +%% ring2z +% Find HEALPix ring number corresponding to a given z coordinate + +%% Syntax +% z = ring2z(nSide,nRing) + +%% Input Arguments +% nSide HEALPix resolution parameter +% nRing numeric array of ring numbers (in range [1,4*nSide-1]) + +%% Return Arguments +% z size(nRing) numeric array of pixel center z coordinates + +%% Example +nRing = reshape(1:15,3,5) +z = ring2z(4,nRing) + +%% See also +% z2ring + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/template_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/template_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/template_help.m (revision 25367) @@ -0,0 +1,64 @@ +%% NAME +% short descr + +%% Syntax + +%% Input Arguments + +%% Return Arguments + +%% Description + +%% Algorithm Notes + +%% Example + +%% See also + +%% Requires + +%% References +% K. M. Gorski et al. 2005. . The Astrophysical Journal 622 759 +% +% K. M. Gorski et al. 1999. . arXiv:astro-ph/9905275 +% +% K. M. Gorski et al. 2002. +% . +% In Astronomical Data Analysis Software and Systems XI, Bohlender, D. A. +% et al. ed. ASP Conference Series Vol 281, pp. 107-111 +% +% Mark R. Calabrettal and Boudewijn F. Roukema. 2007. +% . Monthly Notices of the Royal Astronomical Society 381:2 +% pp. 865 - 872. +% +% W. O'Mullane et al. 2001. . In Mining The Sky. Anthony Banday +% et al., Ed. Anthony ESO Astrophysics Symposia 2001, 638 - 648. +% +% L. Nicastro and Calderone, G. 2008. +% . In Astronomical Data +% Analysis Software and Systems (ADASS) XVII. A. W. Argle et. al, ed. +% Astronomical Society of the Pacific Conference Series vol. 394 pp. +% 486-490. + +%% Notices +% Copyright 2010-2011 Lee Samuel Finn +% Licensed under the Apache License, Version 2.0 (the "License"); +% you may not use this file except in compliance with the License. +% You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, software +% distributed under the License is distributed on an "AS IS" BASIS, +% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +% See the License for the specific language governing permissions and +% limitations under the License. + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/vec2ang_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/vec2ang_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/vec2ang_help.m (revision 25367) @@ -0,0 +1,42 @@ +%% vec2ang +% Convert from cartesian direction vector to position angle on sphere + +%% Syntax +% tp = vec2ang(xyz) + +%% Input Arguments +% xyz cell array of [3,1] cartesian coordinate direction vectors +% +% The direction vectors need not be normalized. + +%% Return Arguments +% tp cell array of [2,1] position angles ([theta, phi]) on sphere + +%% Example +% Create a [3,4] cell array of position angles and manually convert to +% cartesian direction vectors. Use vec2ang to convert back to position +% angles and compare with originals + +% angular position vectors +z = 2*rand(3,4)-1; +th = acos(z); +ph = 2*pi*rand(3,4)-pi; + +% convert to cartesian coordinate vectors +x = sqrt(1-z.^2).*cos(ph); +y = sqrt(1-z.^2).*sin(ph); +xyz = arrayfun(@(x,y,z)([x,y,z]),x,y,z,'UniformOutput',false); + +% use vec2ang to convert cartesian coordinate vectors to angular position +% vectors +tp = vec2ang(xyz); + +% Compare original vectors with final vectors +tp0 = arrayfun(@(t,p)([t;p]),th,ph,'UniformOutput',false); +tf = cellfun(@(x,y)(any(abs(x-y)>4*eps)),tp,tp0) + +%% See also +% ang2vec + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/vec2pix_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/vec2pix_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/vec2pix_help.m (revision 25367) @@ -0,0 +1,37 @@ +%% vec2pix +% Convert cartesian direction vectors to MEALPix pixel numbers + +%% Syntax +% nPix = vec2pix(nSide,xyz,'Param1',Value1,...); + +%% Input Arguments +% nSide HEALPix resolution parameter +% xyz cell array of [3,1] cartesian direction vectors +% +% Param Value +% 'nest' use nested indexing (true | {false}) + +%% Return Arguments +% nPix size(xyz) pixel number array + +%% Example +% Find pixels corresponding to a [4,5] set of (unnormalized) cartesian +% vectors: +nSide = 4; +xyz0 = squeeze(num2cell(2*rand(3,4,5)-1,1)); +nPix = vec2pix(nSide,xyz0) + +% Find the normalized vectors corresponding to the pixels +xyz1 = pix2vec(nSide,nPix); + +% Find angular distance from vectors to pixels +d = angDist(xyz0,xyz1) + +%% See also +% ang2pix, pix2ang, pix2vec + +%% Requires +% ring2nest + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/xy2pix_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/xy2pix_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/xy2pix_help.m (revision 25367) @@ -0,0 +1,37 @@ +%% xy2pix +% Find (nest indexed) MEALPix pixel number corresponding to a pixel index +% on a HEALPix face + +%% Syntax +% ipix = xy2pix(nSide, ix, iy, nFace) + +%% Input Arguments +% nside HEALPix resolution parameter (power of 2) +% ix,iy integer pixel location on face (1 <= ix, iy <= nSide) +% nFace HEALPix face number (1 <= nFace <= 12) +% +% One or both ix and iy may be scalars. If neither is a scalar both must +% have the same size. + +%% Return Arguments +% ipix NESTED pixel index for given x,y on the face + +%% Description +% At lowest resolution HEALPix divides the sphere into 12 faces, each with +% its own coordinate system. Each face is further subdivided in a regular +% fashion as the resolution parameter is increased. Given a face number +% (from 1 to 12), a resolution parameter, and a location on the face XY2PIX +% identifies the pixel number (in the nested indexing scheme). + +%% Example +nSide = 4; +nFace = 3; +ix = [1 2; 1 2]; +iy = [1 1; 2 2]; +nPix = xy2pix(nSide,ix,iy,nFace) + +%% See also +% pix2xy + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/ylm_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/ylm_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/ylm_help.m (revision 25367) @@ -0,0 +1,28 @@ +%% ylm +% Evaluate spherical harmonic basis functions $Y_{LM}$ at pixel centers + +%% Syntax +% v = ylm(nSide,L,M,'Param1',Value1,'Param2',Value2,...); + +%% Input Arguments +% nSide HEALPix resolution parameter (power of 2) +% L spherical harmonic degree (0 <= L) +% M spherical harmonic order (-M <= L <= M) +% +% Param Value +% 'real' return real values (true | {false}) +% 'nest' nested indexing flag (true | {false}) + +%% Return Arguments +% v spherical harmonics evaluated at pixel centers + +%% Example +% Plot Ylm for (L,M) = (4,2) on nSide = 16 HEALPix sphere +pix = ylm(16,4,2); +hp3d(pix); + +%% Requires +% pix2ang + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/help/z2ring_help.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/help/z2ring_help.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/help/z2ring_help.m (revision 25367) @@ -0,0 +1,27 @@ +%% z2ring +% Return HEALPix ring number for a given z + +%% Syntax +% zNum = z2ring(nSide,z) + +%% Input Arguments +% nSide HEALPix resolution parameter +% z numeric array of z values + +%% Return Arguments +% nRing size(z) numeric array of ring numbers + +%% Description +% HEALPix pixels centers are arranged in rings of constant z. Each ring +% spans a given range of z. z2ring returns the ring numbers that cover the +% input z values. + +%% Example +nSide = 4; +zNum = z2ring(nSide,reshape((-7:7)*2/15,3,5)) + +%% See also +% ring2z + +%% +% Copyright 2010-2011 Lee Samuel Finn. . Index: /issm/trunk-jpl/externalpackages/mealpix/info.xml =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/info.xml (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/info.xml (revision 25367) @@ -0,0 +1,61 @@ + + + + + + + + + 2010b + + MEALPix + + + toolbox + + + + ./help + + $toolbox/matlab/icons/bookicon.gif + + + + + + + + + + + + web ./helpfiles/mytoolbox_product_page.html -helpbrowser + + $toolbox/matlab/icons/bookicon.gif + + + + + + + + + + + + + + + + + + + + + web http://gwastro.org -browser; + $docroot/techdoc/matlab_env/examples/webicon.gif + + + + Index: /issm/trunk-jpl/externalpackages/mealpix/src/Contents.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/Contents.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/Contents.m (revision 25367) @@ -0,0 +1,28 @@ +% MEALPix +% alm2pix - Evaluate spherical harmonic function expansion on HEALPix pixels +% alm2spec - valuate angular power spectral density +% ang2pix - Convert spherical coordinate location to HEALPix pixel number +% ang2vec - Convert from spherical to cartesian coordinates +% angDist - Computes angular distance on the unit sphere +% corners - Find pixel corners +% hp3d - Plots a HEALPix dataset on a 3D sphere +% hpInterp - Interpolate data on sphere to MEALPix pixels +% hpStat - Compute descriptive statistics of a pixel value list +% inRing - return pixels in ring or ring section +% nPix2nSide - Convert number of pixels to number of facet side subdivisions +% nSide2nPix - Convert resolution parameter to number of pixels +% nSide2res - Calculates angular resolution of a HEALPix map in arcsec +% nest2ring - Convert MEALPix pixel numbers from nest to ring indexing +% pix2alm - Find spherical harmonic decomposition of function on sphere +% pix2ang - Find spherical coordinate location of HEALPix pixels +% pix2vec - convert HEALPix pixel numbers to (x,y,z) unit vector +% pix2xy - Find pixel cartesian face coordinates +% queryDisc - Find all pixels whose centers are within a specified disc +% res2nSide - Finds minimum base pixel sub-division for required resolution +% ring2nest - Convert MEALPix pixel numbers from ring to nest indexing +% ring2z - Find HEALPix ring number corresponding to a given z coordinate +% vec2ang - Convert cartesian direction vector to position angle on sphere +% vec2pix - Convert cartesian direction vectors to MEALPix pixel numbers +% xy2pix - gives the pixel number ipix (NESTED) of ix, iy and face_num +% ylm - returns spherical harmonic basis Y_L^M on the HEALPix sphere +% z2ring - Find HEALPix ring numbers from direction vector z coordinate Index: /issm/trunk-jpl/externalpackages/mealpix/src/alm2pix.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/alm2pix.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/alm2pix.m (revision 25367) @@ -0,0 +1,76 @@ +function v = alm2pix(nSide, alm, varargin) +% ALM2PIX Evaluate spherical harmonic function expansion on HEALPix pixels +% +% v = alm2pix( nSide, alm, 'Param1', Value1, 'Param2', Value2, ...) +% +% nSide HEALPix resolution parameter (power of 2) +% alm Spherical harmonic expansion coefficients +% +% Param Value +% 'lmax' max order of harmonic to calculate (default floor2*nSide/3) +% 'nest' (true | {false}) +% +% v array of function values at pixels (numel(v) = 12*nSide^2) +% +% Coefficient of Y_{lm} is alm((l+1)*(l+1)+m-1) +% +% Example +% % Shows the gibbs-like error at the poles +% hp3d(alm2pix(16,pix2alm(16,ones(1,12*16^2)))-1) +% +% Requires ylm +% See also pix2alm +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: alm2pix.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +%% Parse and validate input +p = inputParser; +p.addRequired('nSide',@qNSide); +p.addRequired('alm',@qALM); +p.addParamValue('lMax', 2*floor(nSide/3),@qLMax); %% SA corrected it on 11/03/2014 +%p.addParamValue('lMax', floor(2*nSide/3),@qLMax); +p.addParamValue('nest', false, @islogical) + +p.parse(nSide, alm, varargin{:}); + +nSide = p.Results.nSide; +alm = p.Results.alm; +lMax = p.Results.lMax; +nest = p.Results.nest; + +%% +nPix = nSide2nPix(nSide); +v = zeros(1,nPix); +p = 0; % alm index +lMax = max(lMax,sqrt(numel(alm))-1); +for L = 0:lMax + for M = -L:L + v = v + alm(1+p).*ylm(nSide,L,M,'nest',nest); + p = p + 1; + end +end + +return + +function tf = qLMax(lMax) +% QLMAX validate lMax +% +% lMax is a non-negative integer-valued scalar + +tf = isnumeric(lMax) && isscalar(lMax) && floor(lMax) == lMax && lMax >= 0; + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/alm2spec.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/alm2spec.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/alm2spec.m (revision 25367) @@ -0,0 +1,62 @@ +function cL = alm2spec(varargin) +% ALM2SPEC valuate angular power spectral density +% +% cL = alm2spec(alm) +% cL = alm2spec(alm, lMax) +% +% alm real valued coefficients +% lMax max order of harmonic to calculate (default 20) +% +% cL angular power spectrum +% +% Coefficients alm are arranged such that alm(k) is the coefficient +% corresponding to k = (L+1)^2 + M - L +% +% Example +% % Power spectrum estimate of dummy data +% pix = ylm(16,4,2) + ylm(16,7,2) + ylm(16,16,-4) + .1.*rand(1,12*16*16); +% alm = pix2alm(16,pix,100); +% cl = alm2spec(alm,100); +% plot(cl); +% +% See also pix2alm, alm2pix +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: alm2spec.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +%% Parse and validate input +p = inputParser; +p.addRequired('alm',@qAlm); +p.addOptional('lMax',20,@qLMax); +p.parse(varargin{:}); + +alm = p.Results.alm; +lMax = p.Results.lMax; + +%% Compute power spectrum +% Angular power spectrum for given L is sum over m of (|alm|^2/(2L+1)) + +% Preallocate +cL=zeros(1,lMax+1); + +cL(1) = abs(alm(1))^2; +for L = 1:lMax + % Coefficients for L run from L^2+1 to (L+1)^2 + cL(L+1) = sum(abs(alm(L^2+(1:(2*L+1)))).^2); +end +cL = cL./(2*(0:L)+1); + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/ang2pix.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/ang2pix.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/ang2pix.m (revision 25367) @@ -0,0 +1,175 @@ +function pix = ang2pix(varargin) +% ANG2PIX Convert spherical coordinate location to HEALPix pixel number +% +% pix = ang2pix(nSide, v); +% pix = ang2pix(nSide, v, nest); +% +% nSide HEALPix resolution parameter +% v (theta, phi) pairs in radians as either size(v,1) == 2 or as +% cellarray of [2,1] +% nest (optional) nested pixel number scheme (true | {false}) +% +% pix array of pix numbers: +% v on input output +% iscell(v) size(v) cell array of pixel numbers +% isnumeric(v) numeric array of pixel numbers +% +% Example +% vTheta = acos(2*rand(3,5)-1); +% vPhi = 2*pi*rand(3,5); +% v = arrayfun(@(x,y)([x,y]),vTheta,vPhi,'UniformOutput',false); +% nSide = 8; +% pix = ang2pix(nSide,v); +% +% NB: Matlab indexes pixels from 1 +% +% See also vec2pix, pix2ang, pix2vec +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: ang2pix.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +%% Parse and validate arguments +p = inputParser; +p.addRequired('nSide',@qNSide); +p.addRequired('v',@qV); +p.addOptional('nest',false,@islogical) +p.parse(varargin{:}); +nSide = p.Results.nSide; +v = p.Results.v; +nest = p.Results.nest; + +%% unpack v +if isnumeric(v) + ns = size(v); + ns = ns(2:end); + vt = reshape(v(1,:),ns); + vp = reshape(v(2,:),ns); +else + vt = cellfun(@(x)(x(1)),v); + vp = cellfun(@(x)(x(2)),v); +end + +%% Compute pix numbers +pix = f_ang2pix_ring(nSide,vt,vp); + +%% If required convert to nest indexing +if nest + pix = ring2nest(nSide,pix); +end + +return + +%% qV +function tf = qV(v) +% QV test v for inputParser + +if isnumeric(v) + tf = 2 == size(v,1); +elseif iscell(v) + nSize = cellfun('prodofsize',v); + tf = all(2 == nSize(:)); +else + tf = false; +end + +return + +%% f_ang2pix_ring +function ipix = f_ang2pix_ring(nside, theta, phi) + +%======================================================================= +% renders the pixel number ipix (RING scheme) for a pixel which contains +% a point on a sphere at coordinates theta and phi, given the map +% resolution parameter nside +%======================================================================= + +%% f90 I/O Check/Init +if(nside<1 || nside>8192) + msgid = [mfilename ':nSide']; + error(msgid,'%s nSide (%d) out of range',msgid,nside); +end; +if(theta<0 | theta>pi) %#ok + msgid = [mfilename ':theta']; + error(msgid,'%s %e is out of range [0,pi)',msgid,theta); +end; + +z = cos(theta); +za = abs(z); +% in [0,4) +tt = mod(phi,2*pi)/(pi/2); + +ip = zeros(size(z)); +ipix = zeros(size(z)); + +% Logical arrays of pix regions (to remove if statement and vectorize calc) +emask = za <= 2/3; +pmask = ~emask; + +%% Equatorial Region Calculation +temp1(emask) = nside*(0.5+tt(emask)); +temp2(emask) = nside*0.75*z(emask); + +% index of ascending edge line +jp(emask) = fix(temp1(emask)-temp2(emask)); + +% index of descending edge line +jm(emask) = fix(temp1(emask)+temp2(emask)); + +% in {1,2n+1} (ring number counted from z=2/3) +ir = zeros(size(emask)); +ir(emask) = fix(nside + 1 + jp(emask) - jm(emask)); + +% kshift=1 if ir even, 0 otherwise +kshift(emask) = fix(1 - mod(ir(emask),2)); + +%% What is this all about? +nl4 = 4*nside; +% in {0,4n-1} +ip(emask) = fix(fix((jp(emask)+jm(emask)-nside+kshift(emask)+1)/2)); +% if(ip(q) >= nl4) +ip(emask & (ip >= nl4)) = fix(ip(emask & (ip >= nl4)) - nl4); +% end; +ipix(emask) = fix(2*nside*(nside-1) + nl4*(ir(emask)-1) + ip(emask)); + +%% NorthSouth Polar Cap Calculation +%else; +%tt=mod(tt,1); +tp(pmask) = tt(pmask) - fix(tt(pmask)); +tmp(pmask) = nside*sqrt(3*(1-za(pmask))); +% increasing edge line index +jp(pmask) = fix(tp(pmask).*tmp(pmask)); +% decreasing edge line index +jm(pmask) = fix((1-tp(pmask)).*tmp(pmask)); +% ring number counted from the closest pole +ir = zeros(size(pmask)); +ir(pmask) = fix(jp(pmask) + jm(pmask) + 1); +% in {0,4*ir-1} +ip(pmask) = fix(fix(tt(pmask).*ir(pmask))); +newmask1 = (pmask&(ip >= 4*ir)); +%if(ip(pmask) >= 4.*ir(pmask)) +ip(newmask1) = fix(ip(newmask1) - 4*ir(newmask1)); +%end; +newmask2 = (pmask&(z>0)); +%if(z(pmask)>0.) +ipix(newmask2) = fix(2*ir(newmask2).*(ir(newmask2)-1) + ip(newmask2)); +newmask3 = (pmask&~(z>0)); +%else; +nPix = 12*nside^2; +ipix(newmask3) = fix(nPix-2*ir(newmask3).*(ir(newmask3)+1)+ip(newmask3)); +%end; + +ipix = ipix+1; %MEALPix numbering +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/ang2vec.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/ang2vec.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/ang2vec.m (revision 25367) @@ -0,0 +1,61 @@ +function xyz = ang2vec(tp) +% ANG2VEC Convert from spherical to cartesian coordinates +% +% xyz = ang2vec(tp) +% +% tp cell of [2,1] angular locations [theta; phi] in radians +% +% xyz size(tp) celll array of [3,1] vectors +% +% Example +% tp = [acos(2*rand(1,12)-1); 2*pi*rand(1,12)]; +% tp = mat2cell(tp,2,ones(12,1)); +% tp = reshape(tp,3,4); +% xyz = ang2vec(tp); +% xyz{2,3} +% tp{2,3} +% +% See also vec2ang +% +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: ang2vec.m 5831 2011-02-25 23:55:19Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + + +%% IO check +error(nargchk(1,1,nargin,'struct')); +error(nargoutchk(1,1,nargout,'struct')); + +xyz = cellfun(@dAng2Vec,tp,'UniformOutput',false); + +return + +function xyz = dAng2Vec(tp) +% dAng2Vec location on sphere to unit cartesian direction vector +% +% xyz = dAng2Vec(tp) +% +% tp [theta, phi] in radians +% +% xyz [3,1] + +t = tp(1); +z = cos(t); +s = sin(t); +p = tp(2); +xyz = [cos(p)*s;sin(p)*s;z]; + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/angDist.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/angDist.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/angDist.m (revision 25367) @@ -0,0 +1,179 @@ +function rad = angDist(varargin) +%ANGDIST Computes angular distance on the unit sphere +% +% rad = angDist(v0,v1,'Param1',Value1,'Param2',Value2,...) +% +% v0, v1 (numerical or cell) array of points on sphere +% +% Param Value +% 'nSide' HEALPix resolution parameter (integer power 2) +% 'nest' nested indexing (true | {false}) +% +% rad array of angular distances (radians) +% +% Description +% angDist computes the angular distances between points on the sphere. +% Points may be specified as MEALPix pixel numbers, spherical coordinate +% location vectors, or cartesian coordinate vectors. Cartesian coordinate +% input vectors need not be normalized. +% +% Either or both of v0, v1 may be a scalar; otherwise they must both be +% the same size. When one is a scalar than the distance is computed between +% that point and all the points specified in the complementary argument. +% When pix0 and pix1 are arrays then the distances are returned between the +% corresponding points in the two arrays. +% +% If any location is specified as a pixel number then nSide must be +% specified. The indexing scheme for pixels center locations defaults to +% ring; to specify nested indexing set the 'nest' parameter to true. +% +% v0 or v1 specified as numeric arrays: +% If all points in v0 are pixel numbers then v0 may be specified as a +% numeric array (similarly v1, mutatis mundis). nSide must be specified. +% Pixels numbers may be specified in ring indexing or nested indexing. For +% nested indexing the nest parameter must be set to true. +% +% pix0 or pix1 specified as cell arrays: +% Either or both pix0 or pix1 may be specified as a cell array, in which +% case each cell specifies a location on the sphere as either a pixel +% number, a two-component ([theta; phi]) spherical coordinate vector, or a +% three component ([x;y;z]) cartesian vector. nSide must be specified if +% any location is identified as a pixel center. +% +% Example +% % angular distance from pole to equator: +% v0 = {[0;0;1]}; % north pole +% v1 = {[1;0;0]}; % point on equator +% rad = angDist(v0, v1); +% +% % angular distance from pole to pixels 1, 17, 33 for nSide = 4: +% v0 = {[0;0;1]}; +% rad = angDist([1, 17, 33], v0, 'nSide', 4); +% +% % angular distance from pixels 1 to 4, 2 to 5, and 3 to 6 for nSide = 8: +% rad = angDist(8, [1, 2, 3], [4, 5, 6]); +% +% Requires pix2vec, ang2vec +% +% See also queryDisc +% +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: angDist.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + + +%% Verify and parse arguments +error(nargoutchk(1,1,nargout,'struct')); + +p = inputParser; +p.addRequired('v0',@(x)(iscell(x) || isnumeric(x))); +p.addRequired('v1',@(x)(iscell(x) || isnumeric(x))); +p.addParamValue('nest',false,@islogical); +% nSide scalar non-negative integer +p.addParamValue('nSide',0,... + @(x)(isscalar(x) && isnumeric(x) && floor(abs(x)) == x)); +p.parse(varargin{:}); + +v0 = p.Results.v0; +v1 = p.Results.v1; +nest = p.Results.nest; +nSide = p.Results.nSide; + +if ~(1 == numel(v0) || 1 == numel(v1)) + nV0 = size(v0); + nV1 = size(v1); + assert(all(nV0 == nV1),[mfilename ':vDim']); +end + +%% +if isnumeric(v0) + v0 = pix2vec(nSide,v0); +else + v0 = cellfun(@(x)(v2vec(x,nSide,nest)),v0,'UniformOutput',false); +end +if isnumeric(v1) + v1 = pix2vec(nSide,v1); +else + v1 = cellfun(@(x)(v2vec(x,nSide,nest)),v1,'UniformOutput',false); +end + +if 1 == numel(v0) + rad = cellfun(@(x)(f_angdist(v0{1},x)),v1); +elseif 1 == numel(v1) + rad = cellfun(@(x)(f_angdist(x,v1{1})),v0); +else + rad = cellfun(@(x,y)(f_angdist(x,y)),v0,v1); +end + +return + +function dist = f_angdist(v1, v2) +% use pix_tools +%======================================================================= +% call angdist(v1, v2, dist) +% computes the angular distance dist (in rad) between 2 vectors v1 and v2 +% in general dist = acos ( v1 . v2 ) +% except if the 2 vectors are almost aligned. +%======================================================================= +% persistent diff r1 r2 sprod vdiff ; + +% if isempty(r1), r1=zeros(1,3); end; +% if isempty(r2), r2=zeros(1,3); end; +% if isempty(vdiff), vdiff=zeros(1,3); end; +% if isempty(diff), diff=0; end; +% if isempty(sprod), sprod=0; end; +%======================================================================= +% normalize both vectors +r1 = v1/norm(v1); +r2 = v2/norm(v2); +sprod = dot(r1, r2); +if(sprod > 0.999) + % almost colinear vectors + vdiff = r1 - r2; + % norm of difference + diff = sqrt(dot(vdiff,vdiff)); + dist = 2*asin(diff/2); +elseif(sprod < -0.999) ; + % almost anti-colinear vectors + vdiff = r1 + r2; + % norm of sum + diff = sqrt(dot(vdiff,vdiff)); + dist = pi - 2*asin(diff/2); +else + % other cases + dist = acos( sprod ); +end; +return + +function v = v2vec(v,nSide,nest) +% V2VEC pixel, angular locations or vectors to vectors + +switch numel(v) + case 1, % pixel + v = pix2vec(nSide,v,'nest',nest); + case 2, % angular coordinate location + z = cos(v(1)); + s = sin(v(1)); + v = [s*cos(v(2));s*sin(v(2));z]; + case 3, + % nop + otherwise, + msgid = [mfilename ':v2vec']; + error(msgid,msgid); +end +v = reshape(v,[],1); + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/corners.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/corners.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/corners.m (revision 25367) @@ -0,0 +1,75 @@ +function xyz = corners(nSide,varargin) +% CORNERS Find pixel corners +% +% xyz = corners(nSide,nPix,'Param1',Value1,...); +% +% nSide HEALPix resolution parameter +% nPix (optional) pixel list. Default all pixels. +% +% Param Value +% 'nest' nested indexing (true | {false}) +% +% xyz size(nPix) cell array of [3,1] cartesian vector pixel corner +% locations +% +% Example: +% nSide = 2^2; +% % Corners for all 48 nSide=4 pixels in ring-indexed order +% xyzR = corners(nSide); +% size(xyzR) +% xyzR{1} +% xyzR{end-1} +% +% nPix = randi(nSide2nPix(nSide), 3, 4); +% % Corners for an assortment of pixels in nested index order +% xyzN = corners(nSide,nPix,'nest',true); +% % Convert to pixel list to ring indexing +% r = xyzR{nest2ring(nSide,nPix(2,3))} +% % Corners for same pixels requested via ring indexing +% n = xyzN{2,3} +% % Compare +% all(r == n) +% +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: corners.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +%% Validate and parse input +p = inputParser; +p.addRequired('nSide',@qNSide); +p.addOptional('nPix',[],@isnumeric); +p.addParamValue('nest',false,@islogical); + +p.parse(nSide,varargin{:}); +nPix = p.Results.nPix; +nest = p.Results.nest; + +%% Get corners +if isempty(nPix) + nPix = 1:nSide2nPix(nSide); +end +if nest + % convert to ring indexing + nPix = nest2ring(nSide,nPix); +end +nPix = nPix - 1; % convert from MEALPix to HEALPix indexing +[~,xyzc]=f_pix2vec_ring(nSide,nPix(:)); + +xyz = reshape(xyzc,12,[]); +xyz = reshape(num2cell(xyz,1),size(nPix)); +xyz = cellfun(@(x)(reshape(x,3,4)),xyz,'UniformOutput',false); + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/hp3d.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/hp3d.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/hp3d.m (revision 25367) @@ -0,0 +1,88 @@ +function varargout = hp3d(vPix, varargin) +% HP3D Plots a HEALPix dataset on a 3D sphere +% +% h = hp3d(vPix,'Param1',Value1,'Param2',Value2,...); +% +% vPix vector of values at HEALPix pixels +% +% Param Value +% 'nest' nest indexing (true | {false}) +% +% h (optional) array of patch handles for the sphere +% +% Example +% % plot nested pixel numbers on a sphere +% hp3d(1:768,'nest',true); +% +% Requires corners, nPix2nSide +% +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: hp3d.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +%% Parse and validate input +p = inputParser; +p.addRequired('vPix',@qVPix); +p.addParamValue('nest',false,@(x)(isscalar(x) && islogical(x))); + +p.parse(vPix,varargin{:}); +nest = p.Results.nest; + +nSide = nPix2nSide(numel(vPix)); +vPix = reshape(vPix,[],1); + +%% Get pixel corners +c = corners(nSide,'nest',nest); + +%% Plot patches +patchPlot = @(c,v)(... + patch(c(1,:),c(2,:),c(3,:),v,'EdgeColor','none','FaceColor','flat')); + % SA :: patch(c(1,:),c(2,:),c(3,:),v,'EdgeColor','black','FaceColor','white')); +h = cellfun(patchPlot,c,num2cell(reshape(vPix,1,[]))); + +% Cleanup the plot +%view(3); +view([90 90]); +axis equal off +colormap(jet(1024)); +colorbar('southoutside'); + +% slower than opengl, but it doesnt crash X... +set(gcf, 'renderer', 'painters') + +if 1 == nargout + varargout{1} = h; +end + +return + +function tf = qVPix(vPix) +% QVPIX validate vPix +% +% vPix is a numeric vector of length 12*k^2 for k a non-negative integer +% power of 2 + +tf = isnumeric(vPix) && max(size(vPix)) == numel(vPix); + +if tf + nSide = floor(sqrt(numel(vPix)/12)); + ln2 = log2(nSide); + tf = tf && ... + numel(vPix) == 12*nSide^2 && ... + floor(ln2) == ln2; +end + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/hpStat.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/hpStat.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/hpStat.m (revision 25367) @@ -0,0 +1,56 @@ +function stats = hpStat(v) +% HPSTAT Compute descriptive statistics of a pixel value list +% +% stats = hpStat(v) +% +% v array of HEALPix pixel values +% +% stats struct +% field contents +% ntot total number of data points +% nvalid number n of valid data points +% mind minimum of valid data +% maxd maximum of valid data +% average average of valid points +% absdev absolute deviation +% var variance +% stddev standard deviation +% skew skewness factor +% kurt kurtosis factor +% +% Example +% % display stats of Y_3^1 +% v = ylm(32,3,1,'real',true); +% stats = hpStat(abs(v)) +% +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: hpStat.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +stats.ntot = numel(v); +nValid = sum(~isnan(v)); +stats.nValid = nValid; +v = v(~isnan(v)); +stats.mind = min(v); +stats.maxd = max(v); +stats.mean = mean(v); +stats.absdev = mean(abs(v-mean(v))); +stats.var = var(v); +stats.stddev = std(v); +stats.skew = skewness(v); +stats.kurt = kurtosis(v); + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/inRing.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/inRing.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/inRing.m (revision 25367) @@ -0,0 +1,229 @@ +function pixList = inRing(nSide, nRing, varargin) +% INRING return pixels in ring or ring section +% +% pixList = inRing(nSide, nRing, phi0, dPhi, 'Param1', Value1, ...); +% +% nSide HEALPix resolution parameter +% nRing Ring number (1 <= nRing <= 4*nSide-1) +% phi0, dPhi (optional) Ring section longitudes [phi-dPhi,phi+dPhi] +% +% Param Value +% 'nest' return pixels with nested indexing (true | {false}) +% 'strict' return only pixels whose center is in longitude range +% (true | {false}) +% +% pixList size(nRing) cell array of MEALPix pixels in ring defined by +% nRing, phi0, dPhi +% +% nRing may be a numeric array, in which case phi0, dPhi and nest may each +% be either scalar or a size(nRing) array. If strict is false (default) +% then all pixels intersected by longitude range are included in pixList. +% +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: inRing.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + + +%% Parse and validate arguments +p = inputParser; +p.addRequired('nSide',@qnSide); +p.addRequired('nRing',@(n)(qNRing(nSide,n))); +p.addOptional('phi0',0,@isnumeric); +p.addOptional('dPhi',2*pi,@isnumeric); +p.addParamValue('nest',false,@islogical); +p.addParamValue('strict',false,@islogical); + +p.parse(nSide,nRing,varargin{:}); + +phi0 = p.Results.phi0; +dPhi = p.Results.dPhi; +nest = p.Results.nest; +strict = p.Results.strict; + +% Validate dimensions +szRing = size(nRing); + +if 1 == prod(szRing) + if all(1 ~= [numel(phi0), numel(dPhi), numel(nest), numel(strict)]) + msgid = [mfilename ':dims']; + error(msgid,msgid); + end +else + if 1 < numel(phi0) + assert(all(size(phi0) == szRing),[mfilename ':phi0Dims']); + else + phi0 = repmat(phi0,szRing); + end + if 1 < numel(dPhi) + assert(all(size(dPhi) == szRing),[mfilename ':dPhiDims']); + else + dPhi = repmat(dPhi,szRing); + end + if 1 < numel(nest) + assert(all(size(nest) == szRing),[mfilename ':nestDims']); + else + nest = repmat(nest,szRing); + end + if 1 < numel(strict) + assert(all(size(strict) == szRing),[mfilename ':strictDims']); + else + strict = repmat(strict,szRing); + end +end + +%% Initializations +npix = nSide2nPix(nSide); +ncap = 2*nSide*(nSide-1); +pixList = cell(size(nRing)); + +%% Put phi0, dPhi in canonical form +dPhi = abs(dPhi); +phi_low = mod(phi0 - dPhi,2*pi); +phi_hi = mod(phi0 + dPhi,2*pi); +mLong = dPhi >= pi; +if 1 == numel(mLong) + mLong = repmat(mLong,size(nRing)); +end + +%% Identify ring number +nr = zeros(size(nRing)); % ring number +iPix1 = zeros(size(nRing)); % lowest pixel number in the ring +iPix2 = zeros(size(nRing)); % highest pixel number in the ring +kShift = zeros(size(nRing)); + +%%% equatorial region rings +mask = nSide <= nRing & nRing <= 3*nSide; +mir = nRing(mask) - nSide + 1; +iPix1(mask) = ncap + 4*nSide*(mir-1); +iPix2(mask) = iPix1(mask) + 4*nSide - 1; +kShift(mask) = mod(mir,2); +nr(mask) = 4*nSide; + +%%% north pole rings +mask = nRing < nSide; +mir = nRing(mask); +iPix1(mask) = 2*mir.*(mir-1); +iPix2(mask) = iPix1(mask) + 4*mir - 1; +kShift(mask) = 1; +nr(mask) = 4*mir; + +%%% south pole rings +mask = nRing > 3*nSide; +mir = 4*nSide - nRing(mask); +iPix1(mask) = npix - 2*mir.*(mir+1); +iPix2(mask) = iPix1(mask) + 4*mir - 1; +kShift(mask) = 1; +nr(mask) = 4*mir; + +%% Construct pixel list: Full ring case +if any(mLong(:)) + %%% Save all ring pixels: ring indexing + temp = arrayfun(... + @(i1,i2)(i1:i2),... + iPix1(mLong),... + iPix2(mLong),... + 'UniformOutput',false); + [pixList{mLong}] = deal(temp{:}); +end + +%% Construct pixel list: Partial ring case +if any(~mLong(:)) + shift = kShift/2; + + nip = sum(~mLong(:)); + ip_low = zeros(nip,1); + ip_hi = zeros(nip,1); + nList = false(nip,1); + + if any(~strict(:)) + % ~strict = conservative : include every intersected pixel, + % even if pixel CENTER is not in the range [phi_low, phi_hi] + mask = ~strict & ~mLong; + mnr = nr(mask); + mip_low = round(mnr.*phi_low(mask)/(2*pi) - shift(mask)); + mip_hi = round(mnr.*phi_hi(mask)/(2*pi) - shift(mask)); + ip_low(~strict(~mLong)) = mod(mip_low, mnr);% in {0,nr-1} + ip_hi(~strict(~mLong)) = mod(mip_hi, mnr); % in {0,nr-1} + end + + if any(strict) + mask = strict & ~mLong; + mnr = nr(mask); + % strict : include only pixels whose CENTER is in [phi_low, phi_hi] + mip_low = ceil(mnr.*phi_low(mask)/(2*pi) - shift(mask)); + mip_hi = floor(mnr.*phi_hi(mask)/(2*pi) - shift(mask)); + df = mod(mip_low - mip_hi, mnr); + df(df<0) = df(df<0) + mnr(df<0); + + % Mark intervals so small (and away from pixel center) that no pixel + % is included in it + n0 = (df == 1) & (dPhi(mask).*mnr < pi); + nList(n0) = true; + + m0 = mip_low >= mnr; + mip_low(m0) = mip_low(m0) - mnr(m0); + m0 = mip_hi < 0; + mip_hi(m0) = mip_hi(m0) + mnr(m0); + + ip_low(strict(~mLong)) = mip_low; + ip_hi(strict(~mLong)) = mip_hi; + end; + + ip_low = ip_low + iPix1(~mLong(:)); + ip_hi = ip_hi + iPix1(~mLong(:)); + + % identify rings with "wrap-around" longitudinal slices + mTop = ip_low > ip_hi; + ip2 = iPix2(~mLong(:)); + ip1 = iPix1(~mLong(:)); + + topList = arrayfun(... + @(pLow,p1,p2,pHi)([pLow:p2, p1:pHi]),... + ip_low(mTop),ip1(mTop),ip2(mTop),ip_hi(mTop), ... + 'UniformOutput',false); + + % rings with simply-connected longitudinal slices + nTopList = arrayfun(... + @(l,h)(l:h),... + ip_low(~mTop),ip_hi(~mTop),... + 'UniformOutput',false); + + mList = cell(sum(~mLong(:)),1); + [mList{mTop}] = deal(topList{:}); + [mList{~mTop}] = deal(nTopList{:}); + [mList{nList}] = deal([]); + + [pixList{~mLong}] = deal(mList{:}); + +end + +%% Convert ring indexed to nest index as required +if any(nest) + [nestList{1:sum(nest(:))}] = deal(pixList{nest}); + nestList = cellfun(... + @(n)(ring2nest(nSide,n)),... + nestList, 'UniformOutput',false); + [pixList{nest}] = deal(nestList{:}); +end + +%% Convert to MEALPix indexing +pixList = cellfun(@(x)(x+1),pixList,'UniformOutput',false); + +return + +function tf = qNRing(nSide,nRing) +tf = all(1 <= nRing(:)) && all(nRing(:) < 4*nSide); +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/nPix2nSide.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/nPix2nSide.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/nPix2nSide.m (revision 25367) @@ -0,0 +1,36 @@ +function nSide = nPix2nSide(nPix) +% NPIX2NSIDE Convert number of pixels to number of facet side subdivisions +% +% nSide = nPix2nSide(nPix); +% +% nPix number of pixels on sphere +% nSide sub-divisions of HEALPix facet side +% +% Example +% nPix = 12*[1:4].^2; +% nSide = nPix2nSide(nPix) +% +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: nPix2nSide.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +nSide = sqrt(nPix/12); +if any(nSide(:) ~= fix(nSide(:))) + msgid = [mfilename ':nPix']; + error(msgid,msgid); +end + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/nSide2nPix.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/nSide2nPix.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/nSide2nPix.m (revision 25367) @@ -0,0 +1,39 @@ +function nPix = nSide2nPix(nSide) +% NSIDE2NPIX Convert resolution parameter to number of pixels +% +% nPix = nSide2nPix(nSide) +% +% nSide integer power-2 sub-divisions of facet sides. May be an array. +% nPix size(nSide) number of pixels on sphere corresponding to nSide. +% +% HEALPix divides the sphere into 12 basic facets. Each facet is further +% divided into nSide sub-divisions for a total of 12*nSide^2 pixels. +% +% Example +% nPix = nSide2nPix([1:4]) +% +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: nSide2nPix.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +lg = log2(nSide(:)); +if ~all(abs(lg) == fix(lg)) + msgid = [mfilename ':nSide']; + error(msgid,msgid); +end +nPix = 12*nSide.^2; + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/nSide2res.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/nSide2res.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/nSide2res.m (revision 25367) @@ -0,0 +1,37 @@ +function res = nSide2res(nSide) +% NSIDE2RES Calculates angular resolution of a HEALPix map in arcsec +% +% angres = nSide2res(nSide) +% +% nSide HEALPix Resolution (or array of nSides) +% +% angres Angular resolution of the given nSide in arcsec +% +% See also res2nSide +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: nSide2res.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + + +if ~(all(nSide > 0) && all(2.^fix(log2(nSide)) == nSide)) + msgid = [mfilename ':nSide']; + error(msgid,msgid); +end +area = 4*pi*(180*3600/pi)^2; +nPix = nSide2nPix(nSide); +res = sqrt(area./nPix); + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/nest2ring.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/nest2ring.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/nest2ring.m (revision 25367) @@ -0,0 +1,57 @@ +function rPix = nest2ring(nSide,nPix) +% NEST2RING Convert MEALPix pixel numbers from nest to ring indexing +% +% rPix = nest2ring(nSide,nPix) +% +% nSide HEALPix resolution parameters +% nPix numeric array of MEALPix nest indexed pixel numbers +% +% rPix size(nPix) numeric array of ring indexed pixel numbers +% +% See also ring2nest +% +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: nest2ring.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + + +%% Parse and validate input +p = inputParser; +p.addRequired('nSide',@qNSide); +p.addRequired('nPix',@(n)(qPix(nSide,n))); +p.parse(nSide,nPix); + +%% Convert +% Save shape +% convert HEALPix indexing as vector +% convert to MEALPix indexing and restore shape + +sz = size(nPix); +rPix = f_nest2ring(nSide,nPix(:)-1); +rPix = reshape(rPix+1,sz); + +return + +function tf = qPix(nSide,nPix) +tf = isnumeric(nPix); +if tf + nPix = nPix(:); + tf = ... + all(fix(abs(nPix)) == nPix) && ... + all(0 < nPix) && ... + all(nSide2nPix(nSide) >= nPix); +end +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/pix2alm.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/pix2alm.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/pix2alm.m (revision 25367) @@ -0,0 +1,67 @@ +function alm = pix2alm(varargin) +% PIX2ALM Find spherical harmonic decomposition of function on sphere +% +% alm = pix2alm(v) +% alm = pix2alm(v, lmax) +% +% v array of pixel values +% lMax (optional) max order of harmonic to calculate +% +% nPix = numel(v), with nPix = 12*nSide^2 for nSide a power of 2. +% lMax defaults to 2*floor(nSide/3) +% +% alm Spherical harmonic expansion coefficients +% +% alm(j) is coefficient corresponding to order L, M with j = (L+1)^2+M-L +% +% Example +% % estimates |a_l^m| of dummy data +% pix = ylm(16,4,2) + ylm(16,7,2) + ylm(16,16,-4) + .1.*rand(1,12*16*16); +% alm = pix2alm(16,pix,100); +% plot(abs(alm)); +% +% Requires ylm +% See also alm2pix +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: pix2alm.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + + +%% Parse and validate input +p = inputParser; +p.addRequired('v',@qPix); +hLMax = @(v)(2*floor(nPix2nSide(numel(v))/3)); +p.addOptional('lMax', hLMax(varargin{1}), @qLMax); + +p.parse(varargin{:}); + +v = p.Results.v; +lMax = p.Results.lMax; + +%% +% Preallocate alm +alm = zeros(1,(lMax+1)^2); + +nSide = sqrt(numel(v)/12); +for L = 0:lMax + ndx = L^2+1; + for M = -L:L + alm(ndx+L+M) = sum(conj(ylm(nSide,L,M)).*v); + end +end +alm = alm*4*pi/numel(v); + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/pix2ang.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/pix2ang.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/pix2ang.m (revision 25367) @@ -0,0 +1,75 @@ +function tp = pix2ang(nSide, varargin) +% PIX2ANG Find spherical coordinate location of HEALPix pixels +% +% tp = pix2ang(nSide,pix,'Param1',Value1,...) +% +% nSide HEALPix resolution parameter (positive integer power of 2) +% pix (optional) array of pixel values to find coordinates of +% +% Param Values +% 'nest' pixels are given in nested indexing (true | {false}) +% +% tp cell array of [theta; phi] pixel locations (radians) +% +% Example +% % Return theta/phi for all pix of nSide=4, nest scheme +% tp = pix2ang(4,'nest') +% +% See also ang2pix, pix2vec, vec2pix +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: pix2ang.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + + +%% I/O Check +p = inputParser; +p.addRequired('nSide', @qNSide); +p.addOptional('pix',-1,@(n)(qPix(nSide,n))); +p.addParamValue('nest',false,@(n)(isscalar(n) && islogical(n))); + +p.parse(nSide,varargin{:}); + +pix = p.Results.pix; +nest = p.Results.nest; + +% Matlab indexes pixels from 1 +if pix < 0 + pix = 1:nSide2nPix(nSide); +end + +%% Work +% Save size +% Convert to HEALPix convention +% Get pixel location angles +% Restore size + +sz = size(pix); +pix = reshape(pix - 1,1,[]); +if nest + [tparray(1,:), tparray(2,:)]=f_pix2ang_nest(nSide,pix); +else + [tparray(1,:), tparray(2,:)]=f_pix2ang_ring(nSide,pix); +end +tp=reshape(num2cell(tparray,1),sz); + +return + +function tf = qPix(nSide,n) +tf = ... + isnumeric(n(:)) && ... + all(n(:) == fix(abs(n(:)))) && ... + all(n(:)<=nSide2nPix(nSide)); +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/pix2vec.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/pix2vec.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/pix2vec.m (revision 25367) @@ -0,0 +1,104 @@ +function varargout = pix2vec(nSide, varargin) +% PIX2VEC convert HEALPix pixel numbers to (x,y,z) unit vector +% +% [x,y,z] = pix2vec(nSide,pix,'Param1',Value1,...) +% xyz = pix2vec(nSide,pix,'Param1',Value1,...) +% +% nSide HEALPix resolution parameter +% pix (optional) numeric array of pixel numbers to find cartesian +% coordinate vectors for +% +% Param Value +% 'nest' indexing scheme is nested (true | {false}) +% +% x,y,z size(pix) numeric array of unit vector coordinate +% components corresponding to pixels pix +% xyz size(pix) cell array of cartesian unit vectors +% corresponding to pixels pix +% +% Example +% % +% nSide = 2^4; +% nPix = nSide2nPix(nSide); +% pix = randi(nPix,2,3); +% xyz = pix2vec(nSide,pix); +% xyzn = pix2vec(nSide,pix,'nest',true); +% +% Requires nest2ring +% Required by angDist, hpInterp, inRing, queryDisc, ringNum +% See also ang2pix, pix2ang, vec2pix +% +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: pix2vec.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + + +%% IO Check +assert(1 == nargout || 3 == nargout,[mfilename ':nargout']); + +%%% Setup parser +p = inputParser; +p.FunctionName = mfilename; +p.StructExpand = true; + +% Argument declarations +p.addRequired('nSide',@qNSide); +p.addOptional('pix',-1,@(n)(isnumeric(n) && all(floor(n(:))==n(:)))); +p.addParamValue('nest',false,@islogical); + +%%% parse and validate input +p.parse(nSide,varargin{:}); +nSide = p.Results.nSide; +pix = p.Results.pix; +nest = p.Results.nest; + +%% Get pixel vector + +if 1 == numel(pix) && pix < 0 + pix = 1:nSide2nPix(nSide); +end +assert(... + all(0 < pix(:)) || all(pix(:) <= nSide2nPix(nSide)),... + [mfilename ':pix']); + +if nest + pix = nest2ring(nSide,pix); +end +xyz = f_pix2vec_ring(nSide,pix-1); + +%% Format return +switch nargout + case 1, + nSize = size(xyz); + if 2 == numel(nSize) + nSize = [nSize 1]; + end + xyz = reshape(xyz,3,[]); + xyz = num2cell(xyz,1); + varargout{1} = reshape(xyz,nSize(2:end)); + case 3, + nSize = size(xyz); + if 2 == numel(nSize); + nSize = [nSize 1]; + end + sz = nSize(2:end); + xyz = reshape(xyz,3,[]); + varargout{1} = reshape(xyz(1,:),sz); + varargout{2} = reshape(xyz(2,:),sz); + varargout{3} = reshape(xyz(3,:),sz); +end + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/pix2xy.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/pix2xy.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/pix2xy.m (revision 25367) @@ -0,0 +1,64 @@ +function xy = pix2xy(nSide, nPix) +%% PIX2XY Find pixel cartesian face coordinates +% +% xy = pix2xy(nSide,nPix) +% +% nSide HEALPix resolution parameter +% nPix nested indexing pixel numbers to find coordinates of +% +% xy size(nPix) cell array of pixel coordinates on face +% +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: pix2xy.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + + +if ~qNSide(nSide) + msgid = [mfilename ':nSide']; + error(msgid,msgid); +end + +if ~(all(nPix(:) > 0) && all(nPix(:) <= nSide2nPix(nSide))) + msgid = [mfilename ':nPix']; + error(msgid,msgid); +end +sPix = size(nPix); +if 2 == numel(sPix) + sPix = [sPix, 1]; +end +ixy = zeros([2,sPix(:)']); + +[pix2x,pix2y]=mk_pix2xy; + +ipf = nPix(:)-1; + +% content of the last 10 bits +ip_low = fix(mod(ipf,1024)); + +% truncation of the last 10 bits +ip_trunc = fix(ipf/1024); + +% content of the next 10 bits +ip_med = fix(mod(ip_trunc,1024)); + +% content of the high weight 10 bits +ip_hi = fix(ip_trunc/1024); +ixy(1,:) = fix(1024*pix2x(ip_hi+1) + 32*pix2x(ip_med+1) + pix2x(ip_low+1)); +ixy(2,:) = fix(1024*pix2y(ip_hi+1) + 32*pix2y(ip_med+1) + pix2y(ip_low+1)); + +xy = reshape(num2cell(ixy,1),sPix); + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/private/correct_ring_phi.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/private/correct_ring_phi.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/private/correct_ring_phi.m (revision 25367) @@ -0,0 +1,27 @@ +function [iring, iphi] = correct_ring_phi(location, iring, iphi) + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +delta = zeros(size(iphi)); +delta(iphi<0) = 1; +delta(iphi>4*iring) = -1; +nzMask = delta ~= 0; +if any(nzMask) + if any(abs(location(nzMask)) ~= 1) + msgid = [mfilename,':location']; + error(msgid,msgid); + end + iring(nzMask) = iring(nzMask) - location(nzMask).*delta(nzMask); + iphi(nzMask) = iphi(nzMask) + delta(nzMask).*(4*iring(nzMask)); +end +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/private/f_nest2ring.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/private/f_nest2ring.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/private/f_nest2ring.m (revision 25367) @@ -0,0 +1,104 @@ +function ipring = f_nest2ring(nSide, ipnest) + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +% $Id: f_nest2ring.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +nPix = nSide2nPix(nSide); +if ~(all(ipnest>=0) && all(ipnest= nSide & (jr <= 3*nSide); +jrE = jr(mask); +nrE = fix(nSide); +n_before(mask) = fix(2*nrE*(2*jrE - nrE - 1)); +kshift(mask) = mod(jrE - nSide, 2); +nr(mask) = nrE; + +% south pole region +mask = jr > 3*nSide; +jrS = jr(mask); +nrS = fix(nl4 - jrS); +n_before(mask) = fix(nPix - 2*nrS.*(nrS + 1)); +kshift(mask) = 0; +nr(mask) = nrS; + +% computes the phi coordinate on the sphere, in [0,2Pi] +% 'phi' number in the ring in {1,4*nr} +jp = fix((jpll(face_num+1).*nr + jpt + 1 + kshift)/2); +maskH = jp>nl4; +maskL = jp<1; +jp(maskH) = fix(jp(maskH) - nl4); +jp(maskL) = fix(jp(maskL) + nl4); + +% in {0, nPix-1} +ipring = fix(n_before + jp - 1); + +return + Index: /issm/trunk-jpl/externalpackages/mealpix/src/private/f_pix2ang_nest.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/private/f_pix2ang_nest.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/private/f_pix2ang_nest.m (revision 25367) @@ -0,0 +1,108 @@ +function [theta, phi] = f_pix2ang_nest(nSide, ipix) + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +% $Id: f_pix2ang_nest.m 6016 2011-04-27 02:52:56Z lsf@GRAVITY.PSU.EDU $ + +% coordinate of the lowest corner of each face +% in unit of nSide +jrll =[ 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4 ]'; +% in unit of nSide/2 +jpll =[ 1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7 ]'; +%----------------------------------------------------------------------- +% total number of points +npix = nSide2nPix(nSide); +if ~(all(ipix>=0) && all(ipix (x,y) mapping +[pix2x,pix2y] = mk_pix2xy; + +npface = nSide^2; +nl4 = 4*nSide; + +% find face and pixel number in the face +% face number in {0,11} +face_num = fix(ipix/npface); +% pixel number in the face {0,npface-1} +ipf = mod(ipix,npface); + +fact1 = 1/(3*nSide^2); +fact2 = 2/(3*nSide); +% finds the x,y on the face (starting from the lowest corner) +% from the pixel number +ix = 0; +iy = 0; +scalemlv = 1; +ismax = 4; +for i=0: ismax; + ip_low = mod(ipf,1024); + ix = fix(ix + scalemlv*pix2x(ip_low+1)); + iy = fix(iy + scalemlv*pix2y(ip_low+1)); + scalemlv = scalemlv*32; + ipf = fix(ipf/1024); +end +ix = fix(ix + scalemlv*pix2x(ipf+1)); +iy = fix(iy + scalemlv*pix2y(ipf+1)); + +%% transforms to (horizontal, vertical) coordinates +% 'vertical' in {0,2*(nSide-1)} +jrt = fix(ix + iy); +% 'horizontal' in {-nSide+1,nSide-1} +jpt = fix(ix - iy); + +%% compute z coordinate on the sphere +% ring number in {1,4*nSide-1} +jr = fix(jrll(face_num+1)*nSide - jrt - 1); + +z = zeros(size(jr)); +nr = zeros(size(jr)); +kshift = zeros(size(jr)); + +% north pole region +mask = jr 3*nSide; +nrM = fix(nl4 - jr(mask)); +z(mask) = - 1 + nrM.^2*fact1; +kshift(mask) = 0; +nr(mask) = nrM; + +%% Convert z to theta +theta = acos(z); + +%% Compute phi in [0,2*pi) +% 'phi' number in the ring in {1,4*nr} +jp =fix(((jpll(face_num+1).*nr) + jpt + 1 + kshift)/2); +if(jp > nl4) + jp = fix(jp - nl4); +end; +if(jp < 1) + jp = fix(jp + nl4); +end; +phi = (pi/2)*(jp -(kshift+1)/2)./nr; + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/private/f_pix2ang_ring.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/private/f_pix2ang_ring.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/private/f_pix2ang_ring.m (revision 25367) @@ -0,0 +1,78 @@ +function [theta, phi] = f_pix2ang_ring(nSide, iPix) +% F_PIX2ANG_RING convert ring-indexed pixel to angular position vector +% +% [theta, phi] = f_pix2ang_ring(nSide,iPix) +% +% nSide HEALPix resolution parameter +% iPix HEALPix ring-indexed pixel number +% +% Hand conversion of HEALPix distribution f90 code pix2ang_ring +% +% Author: Lee Samuel Finn +% Copyring 2011 + +% $Id: f_pix2ang_ring.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + + +npix = nSide2nPix(nSide); +if any(iPix<0) || any(iPix>npix-1) + msgid = [mfilename ':ndx']; + error(msgid,msgid); +end + +nl2 = 2*nSide; + +% points in each polar cap, =0 for nSide =1 +nCap = nl2*(nSide-1); + +% preallocate results +theta = zeros(size(iPix)); +phi = zeros(size(iPix)); + +%% North Polar cap ------------- +mask = iPix < nCap; +% counted from North pole +iPixM = iPix(mask); +iRingM = round(sqrt((iPixM+1)/2)); +iPhiM = fix(iPixM - 2*iRingM.*(iRingM - 1)); +% fix round-off error appearing at large nSide +[iRingM, iPhiM]=correct_ring_phi(1, iRingM, iPhiM); +theta(mask) = acos(1-(iRingM/nSide).^2/3); +phi(mask) = (pi/2)*(iPhiM + 0.5)./iRingM; + +%% Equatorial region ------ +mask = ~mask & (iPix < npix-nCap); +ipM = fix(iPix(mask) - nCap); +nl4 = 4*nSide; +% counted from North pole +iRingM = fix(ipM/nl4) + nSide; +iPhiM = mod(ipM,nl4); +% 0 if iring+nSide is odd, 1/2 otherwise +fodd = 0.5*mod(iRingM+nSide+1,2); +theta(mask) = acos((nl2 - iRingM)/(1.5*nSide)); +phi(mask) = (pi/2)*(iPhiM + fodd)/nSide; + +%% South Polar cap ----------------------------------- +mask = npix-nCap <= iPix; +ipM = fix(npix - iPix(mask)); +% counted from South pole +iRingM = round(sqrt(ipM/2)); +iPhiM = fix(2*iRingM.*(iRingM + 1) - ipM); +% fix round-off error appearing at large nSide +[iRingM, iPhiM]=correct_ring_phi(-1, iRingM, iPhiM); +theta(mask) = acos((iRingM/nSide).^2/3 - 1); +phi(mask) = (pi/2)*(iPhiM + 0.5)./iRingM; + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/private/f_pix2vec_ring.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/private/f_pix2vec_ring.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/private/f_pix2vec_ring.m (revision 25367) @@ -0,0 +1,202 @@ +function varargout = f_pix2vec_ring(nSide, hPix) +% F_PIX2VEC_RING conversion of f90 PIX2VEC_RING HEALPix code +% +% [vector, vertex] = f_pix2vec_ring(nSide,hPix) +% +% nSide (integer) sub-divisions of HEALPix facets +% hPix (vector integers) HEALPix pixel numbers +% (0 <= hPix < 12*nSide^2) +% +% vector cell(size(hPix)) of [3,1] unit vectors to pixel centers +% vertex (optional) cell(size(hPix)) of [3,4] unit vectors to pixel +% vertices (corners) +% +% renders vector (x,y,z) coordinates of the nominal pixel center +% for the pixel number hPix (RING scheme) +% given the map resolution parameter nSide +% also returns the (x,y,z) position of the 4 pixel vertices (=corners) +% in the order N,W,S,E +% +% Hand rewrite of f90 code +% +% Author: Lee Samuel Finn +% Copyright 2010 + +% $Id: f_pix2vec_ring.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +%% I/O Check/Init +error(nargoutchk(1,2,nargout,'struct')); +error(nargchk(2,2,nargin,'struct')); + +nPix = nSide2nPix(nSide); +if ~(... + all(floor(hPix(:)) == abs(hPix(:))) && all(hPix(:) < nPix)) + msgid = [mfilename ':hPix']; + error(msgid,msgid); +end + +% in {1, nPix} +hPix1 = fix(hPix + 1); +nl2 = fix(2*nSide); +nl4 = fix(4*nSide); +% points in each polar cap, = 0 for nSide = 1 +nCap = fix(2*nSide*(nSide-1)); +fact1 = 1.5*nSide; +fact2 = 3.0*nSide^2; + +% Calculate corners (vertices)? +do_vertex = (nargout == 2); + +% Associate pixels with regions (north cap, equatorial, south cap) +ncMask = hPix1 <= nCap; +eqMask = hPix1 <= nl2*(5*nSide+1) & ~ncMask; +scMask = ~(eqMask | ncMask); + +% Pre-allocate +z = zeros(size(hPix)); +phi = zeros(size(hPix)); +hDelta_Phi = zeros(size(hPix)); +z_nv = zeros(size(hPix)); +z_sv = zeros(size(hPix)); +phi_nv = zeros(size(hPix)); +phi_sv = zeros(size(hPix)); + +% North polar cap +if any(ncMask(:)) + hip = hPix1(ncMask)/2; + fihip = fix(hip); + iRing = fix(sqrt(hip-sqrt(fihip)))+1; + iPhi = hPix1(ncMask) - 2*iRing.*(iRing-1); + + z(ncMask) = 1-iRing.^2/fact2; + phi(ncMask) = (iPhi-0.5).*pi./(2*iRing); + + if (do_vertex) + hDelta_Phi(ncMask) = pi./(4*iRing); % half pixel width + z_nv(ncMask) = 1-(iRing-1).^2/fact2; + z_sv(ncMask) = 1-(iRing+1).^2/fact2; + iPhi_Mod = mod(iPhi-1,iRing); % in {0..iRing-1} + iPhi_Rat = fix((iPhi-1)./iRing); % in {0,1,2,3} + phi_nv_m = phi_nv(ncMask); + ndx = (1 < iRing); + phi_nv_m(ndx) = (iPhi_Rat(ndx)+iPhi_Mod(ndx)./(iRing(ndx)-1))*pi/2; + phi_nv(ncMask) = phi_nv_m; + phi_sv(ncMask) = (iPhi_Rat+(iPhi_Mod+1)./(iRing+1))*pi/2; + end +end + +% Equatorial region +if any(eqMask(:)) + ip = hPix1(eqMask) - nCap -1; + iRing = fix(ip/nl4) + nSide; % counted from north pole + iPhi = mod(ip,nl4)+1; + + fOdd = 0.5*(1+mod(iRing+nSide,2)); % 1 (1/2) if iRing+nSide is odd (even) + z(eqMask) = (nl2-iRing)/fact1; + phi(eqMask) = (iPhi-fOdd)*pi/(2*nSide); + + if (do_vertex) + hDelta_Phi(eqMask) = pi./(4*nSide); % half pixel width + phi_nv_m = phi(eqMask); + phi_sv_m = phi(eqMask); + z_nv_m = (nl2-iRing+1)/fact1; + z_sv_m = (nl2-iRing-1)/fact1; + nNdx = (iRing == nSide); % northern transition + sNdx = (iRing == 3*nSide); % southern transition + ndx = sNdx|nNdx; % transition + z_nv_m(nNdx) = 1-(nSide-1)^2/fact2; + z_sv_m(sNdx) = -1+(nSide-1)^2/fact2; + iPhi_Mod(ndx) = mod(iPhi(ndx)-1,nSide); % in {0..nSide-1} + iPhi_Rat(ndx) = fix((iPhi(ndx)-1)/nSide); % in {0..3} + if nSide > 1 + phi_nv_m(nNdx) = (iPhi_Rat(nNdx)+iPhi_Mod(nNdx)/(nSide-1))*pi/2; + phi_sv_m(sNdx) = (iPhi_Rat(sNdx)+iPhi_Mod(sNdx)/(nSide-1))*pi/2; + end + phi_nv(eqMask) = phi_nv_m; + phi_sv(eqMask) = phi_sv_m; + z_nv(eqMask) = z_nv_m; + z_sv(eqMask) = z_sv_m; + end +end + +% South polar cap +if any(scMask(:)) + ip = nPix - hPix1(scMask) + 1; + hip = ip/2; + fihip = fix(hip); + iRing = fix(sqrt(hip-sqrt(fihip)))+1; % counted from south pole + iPhi = 4*iRing+1-(ip-2*iRing.*(iRing-1)); + + z(scMask) = -1+iRing.^2/fact2; + phi(scMask) = (iPhi-0.5).*pi./(2*iRing); + + if (do_vertex) + hDelta_Phi(scMask) = pi./(4*iRing); + z_nv(scMask) = -1+(iRing+1).^2/fact2; + z_sv(scMask) = -1+(iRing-1).^2/fact2; + iPhi_Mod = mod(iPhi-1,iRing); % in {0..iRing-1} + iPhi_Rat = fix((iPhi-1)./iRing); % in {0..3} + phi_nv(scMask) = (iPhi_Rat+(iPhi_Mod+1)./(iRing+1))*pi/2; + ndx = iRing > 1; + phi_sv_m = phi_sv(scMask); + phi_sv_m(ndx) = (iPhi_Rat(ndx)+iPhi_Mod(ndx)./(iRing(ndx)-1))*pi/2; + phi_sv(scMask) = phi_sv_m; + end +end + +sth = sqrt(1-z).*sqrt(1+z); + +%% Vertices +if (do_vertex) + v = zeros(horzcat([3,4],size(hPix))); + + % West vertex + phi_wv = phi - hDelta_Phi; + v(1,2,:) = sth(:).*cos(phi_wv(:)); + v(2,2,:) = sth(:).*sin(phi_wv(:)); + v(3,2,:) = z(:); + + % East vertex + phi_ev = phi + hDelta_Phi; + v(1,4,:) = sth(:).*cos(phi_ev(:)); + v(2,4,:) = sth(:).*sin(phi_ev(:)); + v(3,4,:) = z(:); + + % North vertex + sth_nv = sqrt(1-z_nv).*sqrt(1+z_nv); + v(1,1,:) = sth_nv(:).*cos(phi_nv(:)); + v(2,1,:) = sth_nv(:).*sin(phi_nv(:)); + v(3,1,:) = z_nv(:); + + % South vertex + sth_sv = sqrt(1-z_sv).*sqrt(1+z_sv); + v(1,3,:) = sth_sv(:).*cos(phi_sv(:)); + v(2,3,:) = sth_sv(:).*sin(phi_sv(:)); + v(3,3,:) = z_sv(:); + + varargout{2} = v; +end + +%% Pixel centers +sth = reshape(sth,[],1); +phi = reshape(phi,[],1); +vector = zeros(horzcat(3,size(hPix))); +vector(1,:) = sth.*cos(phi); +vector(2,:) = sth.*sin(phi); +vector(3,:) = z(:); + +varargout{1} = vector; + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/private/f_ring2nest.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/private/f_ring2nest.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/private/f_ring2nest.m (revision 25367) @@ -0,0 +1,135 @@ +function ipnest = f_ring2nest(nside, ipring) + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +% $Id: f_ring2nest.m 6224 2011-08-22 18:49:47Z lsf@GRAVITY.PSU.EDU $ + +npix = nSide2nPix(nside); +assert(all(ipring>=0) && all(ipring mifm); +mface(smask) = fix(mifp(smask)+7); % (half-)faces 8 to 11 + +face_num(mask) = mface; +irn(mask) = mirn; +iphi(mask) = miphi; + +%% south polar cap +mask = ipring >= (npix - ncap); +mRing = ipring(mask); + +mip = fix(npix - mRing); +% counted from South pole +mirs = round(sqrt(mip/2)); +miphi = 2*mirs.*(mirs + 1) - mip; +[mirs, miphi] = correct_ring_phi(1, mirs, miphi); +kshift(mask) = 0; +nr(mask) = fix(mirs); +irn(mask) = fix(nl4 - mirs); +face_num(mask) = fix(miphi./mirs) + 8; +iphi(mask) = miphi; + +%% finds the (x,y) on the face +irt = irn - jrll(face_num+1)*nside + 1; % in {-nside+1,0} +ipt = 2*iphi - jpll(face_num+1).*nr - kshift + 1; % in {-nside+1,nside-1} + +mask = ipt >= nl2; % for the face #4 +ipt(mask) = ipt(mask) - 8*nside; + +ix = fix((ipt - irt )/2); +iy = -fix((ipt + irt )/2); + +scalemlv = 1; +scale_factor = 16384; +ipf = 0; +% for nside in [2^14, 2^20] +ismax = 1; +if (nside > 1048576) + ismax = 3; +end; +for i = 0:ismax + % last 7 bits + ix_low = mod(ix, 128); + iy_low = mod(iy, 128); + + ipf = fix(ipf +(x2pix(ix_low+1)+y2pix(iy_low+1))*scalemlv); + scalemlv = scalemlv*scale_factor; + % truncate out last 7 bits + ix = fix(ix/128); + iy = fix(iy/128); +end +ipf = fix(ipf +(x2pix(ix+1)+y2pix(iy+1))*scalemlv); + +% in {0, 12*nside**2 - 1} +ipnest = fix(ipf + (face_num.*(fix(npix/12)))); +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/private/f_surface_triangle.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/private/f_surface_triangle.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/private/f_surface_triangle.m (revision 25367) @@ -0,0 +1,56 @@ +function surface = f_surface_triangle(vec1, vec2, vec3) + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +% $Id: f_surface_triangle.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% use pix_tools +%======================================================================= +% returns the surface in steradians +% of the spherical triangle with vertices vec1, vec2, vec3 +% +% algorithm : finds triangle sides and uses l'Huilier formula to compute +% 'spherical excess' = surface area of triangle on a sphere of radius one +% see, eg Bronshtein, Semendyayev Eq 2.86 +%======================================================================= +% real(kind=dp), dimension(1:3) :: v1, v2, v3 +% persistent side x0 x1 x2 x3 ; +% +% if isempty(side), side=zeros(1,3); end; +% % real(kind=dp) :: hp +% if isempty(x0), x0=0; end; +% if isempty(x1), x1=0; end; +% if isempty(x2), x2=0; end; +% if isempty(x3), x3=0; end; +%======================================================================= +% half perimeter +% hp = 0.5 * (side1 + side2 + side3) +% ! l'Huilier formula +% x0 = tan( hp * 0.5) +% x1 = tan((hp - side1) * 0.5) +% x2 = tan((hp - side2) * 0.5) +% x3 = tan((hp - side3) * 0.5) +% find triangle sides +side(1)=f_angdist(vec2, vec3); +side(2)=f_angdist(vec3, vec1); +side(3)=f_angdist(vec1, vec2); +% divide by 4 +side([1:3]) = side([1:3]) .* 0.25; +% l'Huilier formula +x0 = tan( side(1) + side(2) + side(3) ); +x1 = tan(-side(1) + side(2) + side(3) ); +x2 = tan( side(1) - side(2) + side(3) ); +x3 = tan( side(1) + side(2) - side(3) ); +surface = 4.0 .* atan( sqrt(x0 .* x1 .* x2 .* x3) ); + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/private/mk_pix2xy.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/private/mk_pix2xy.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/private/mk_pix2xy.m (revision 25367) @@ -0,0 +1,55 @@ +function [pix2x,pix2y] = mk_pix2xy + +%======================================================================= +% constructs the array giving x and y in the face from pixel number for the +% nested (quad-cube like) ordering of pixels +% +% the bits corresponding to x and y are interleaved in the pixel number one +% breaks up the pixel number by even and odd bits +%======================================================================= + +% $Id: mk_pix2xy.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +pix2x=zeros(1024,1); +pix2y=zeros(1024,1); +for kpix=0:1023 + jpix = kpix; + ix = 0; + iy = 0; + % bit position (in x and y) + ip = 1; + while (1); + % go through all the bits + if(jpix == 0) + break + end + % bit value (in kpix), goes in ix + id = mod(jpix,2); + jpix = fix(jpix/2); + ix = id.*ip+ix; + % bit value (in kpix), goes in iy + id = mod(jpix,2); + jpix = fix(jpix/2); + iy = id.*ip+iy; + % next bit (in x and y) + ip = 2*ip; + end + % in 0,31 + pix2x(kpix+1) = ix; + % in 0,31 + pix2y(kpix+1) = iy; +end + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/private/mk_xy2pix.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/private/mk_xy2pix.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/private/mk_xy2pix.m (revision 25367) @@ -0,0 +1,48 @@ +function [x2pix y2pix]=mk_xy2pix +%======================================================================= +% sets the array giving the number of the pixel lying in (x,y) +% x and y are in {1,128} +% the pixel number is in {0,128**2-1} +% +% if i-1 = sum_p=0 b_p * 2^p +% then ix = sum_p=0 b_p * 4^p +% iy = 2*ix +% ix + iy in {0, 128**2 -1} +%======================================================================= + +% $Id: mk_xy2pix.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + + +%for converting x,y into pixel numbers +x2pix=zeros(128,1); +y2pix=zeros(128,1); +for i = 1:128; + j = (i-1); + k = 0; + ip = 1; + while (1); + if(j==0) + x2pix(i) = k; + y2pix(i) = 2*k; + break; + else + id = mod(j,2); + j = fix(j/2); + k = ip*id+k; + ip = 4*ip; + end + end +end +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/private/qALM.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/private/qALM.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/private/qALM.m (revision 25367) @@ -0,0 +1,27 @@ +function tf = qALM(alm) +% QALM validate alm +% +% alm is a numeric vector and numel(alm) = k^2 for some +% postive integer k +% +% Required by alm2spec, alm2pix +% Author: Lee Samuel Finn +% Copyright 2010-2011 + +% $Id: qALM.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +tf = isnumeric(alm) && sqrt(numel(alm))^2 == numel(alm); + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/private/qLMax.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/private/qLMax.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/private/qLMax.m (revision 25367) @@ -0,0 +1,14 @@ +function tf = qLMax(lMax) +% QLMAX validate lMax +% +% lMax is a non-negative integer-valued scalar +% +% Required by alm2spec, alm2pix +% Author: Lee Samuel Finn +% Copyright 2010 + +% $Id: qLMax.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +tf = isnumeric(lMax) && isscalar(lMax) && floor(lMax) == abs(lMax); + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/private/qNSide.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/private/qNSide.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/private/qNSide.m (revision 25367) @@ -0,0 +1,27 @@ +function tf = qNSide(nSide) +% QNSIDE validate nSide +% +% tf = qNSide(nSide) +% +% returns true if nSide is a numeric scalar non-negative integer-valued +% power of 2 + +% $Id: qNSide.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + + +tf = isnumeric(nSide) && isscalar(nSide) && 2^floor(log2(nSide)) == nSide; + +return + Index: /issm/trunk-jpl/externalpackages/mealpix/src/private/qPix.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/private/qPix.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/private/qPix.m (revision 25367) @@ -0,0 +1,23 @@ +function tf = qPix(pix) +% QPIX validate list of pixel values +% +% pix is a numeric vector with length 12*k^2 for integer k + +% $Id: qPix.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +tf = isnumeric(pix) && numel(pix) == max(size(pix)); +tf = tf && floor(sqrt(numel(pix)/12))^2*12 == numel(pix); + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/private/ringNum.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/private/ringNum.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/private/ringNum.m (revision 25367) @@ -0,0 +1,60 @@ +function nRing = ringNum(nSide, z, varargin) +%% RINGNUM Find ring number for given z coordinate +% +% ring = ringNum(nSide, z, shift) +% +% nSide HEALPix resolution parameter (integer power 2) +% z numeric array of z coordinates on unit sphere +% shift (optional) ring number "rounding" mode +% shift = 0: return ring number whose center is closest to z +% shift < 0: return ring number immediately north of closest ring +% shift > 0: return ring number immediately south of closest ring +% +% nRing size(z) numeric array of ring numbers +% +% Author: Lee Samuel Finn +% Copyright: 2011 + +% $Id: ringNum.m 5792 2011-02-20 03:21:25Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + + +%% Parse and validate input +p = inputParser; +p.addRequired('nSide',@qNSide); +p.addRequired('z',@(x)(all(abs(x(:))<=1))); +p.addOptional('shift',0,@isnumeric); +p.parse(nSide,z,varargin{:}); +shift = p.Results.shift/2; + +%% preallocate +nRing = zeros(size(z)); + +%% equatorial regime +mask = abs(z) <= 2/3; +nRing(mask) = round(nSide*(2-1.5*z(mask)) + shift); + +%% north polar cap +mask = z > 2/3; +nTmp = round(nSide*sqrt(3*(1-z(mask))) + shift); +nTmp(0 == nTmp) = 1; +nRing(mask) = nTmp; + +%% south polar cap +mask = z < -2/3; +nTmp = round(nSide*sqrt(3*(1+z(mask))) - shift); % Yes: this is a -shift +nTmp(0 == nTmp) = 1; +nRing(mask) = 4*nSide - nTmp; + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/queryDisc.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/queryDisc.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/queryDisc.m (revision 25367) @@ -0,0 +1,86 @@ +function pix = queryDisc( nSide, c0, r0, varargin ) +% QUERYDISC Find all pixels whose centers are within a specified disc +% +% +% pix = queryDisc(nSide,c0,r,'Param1',Value1,...) +% +% nSide +% c0 disc center. May be specified as a pixel number, anuglar +% position vector [theta, phi] in rads, or as a cartesian unit +% vector [x;y;z]. +% r disc radius (rads) +% +% Param Value +% 'nest' pixel indexing is nested (true | {false}) +% +% pix list of pixels in disc +% +% Example +% nSide=8; +% pix = queryDisc(nSide,1,pi/10,'nest'); +% numel(pix) +% nSide2nPix(nSide)*((pi/10)/(4 Pi)) +% +% See also queryPolygon, queryStrip, queryTriangle +% +% Requires pix2vec, ang2vec, angDist +% +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: queryDisc.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +%% validate and parse arguments +p = inputParser; +p.addRequired('nSide',@qNSide); +p.addRequired('c0',@(x)(isnumeric(x) && any(numel(x) == 1:3))); +p.addRequired('r0',@(x)(isnumeric(x) && isscalar(x) && x > 0 && x <= pi)); +p.addParamValue('nest',false,@(n)(isscalar(n) && islogical(n))); +p.parse(nSide,c0,r0,varargin{:}); +nest = p.Results.nest; + +%% Setup vIn +switch numel(c0) + case 1 %pix + vc = cell2mat(pix2vec(nSide,c0,'nest',nest)); + case 2 %ang + vc = cell2mat(ang2vec(c0)); + case 3 %vec + vc = c0/norm(c0); +end + +%% Computation +% Find max, min z of disk +% Find corresponding ring numbers +% Find pixels in rings +% Select pixels within ring + +th = acos(vc(3)); +zMax = cos(max(th-r0,0)); +zMin = cos(min(th+r0,pi)); + +ringMin = z2ring(nSide,zMax); +ringMax = z2ring(nSide,zMin); + +rings = min(ringMin,ringMax):max(ringMin,ringMax); +pix = arrayfun(@(r)(cell2mat(inRing(nSide,r))),rings,'UniformOutput',false); +pix = cell2mat(pix); + +d = angDist(pix,{vc},'nSide',nSide); +pix = pix(d < r0); % closer than r0 +if nest + pix = ring2nest(nSide,pix); +end +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/res2nSide.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/res2nSide.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/res2nSide.m (revision 25367) @@ -0,0 +1,35 @@ +function nSide = res2nSide(r) +% RES2NSIDE Finds minimum base pixel sub-division for required resolution +% +% nSide = res2nSide(r) +% +% r angular resolution (arc seconds) +% +% nSide minimum base pixel sub-divisions corresponding to resolution r +% +% See also nSide2res +% +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: res2nSide.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + + +nPix0 = 4*pi./(r*pi/(180*3600)).^2; +nSide0 = realsqrt(nPix0/12); +nSide = pow2(ceil(log2(nSide0)-2*eps)); +nSide(nSide<1) = 1; + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/ring2nest.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/ring2nest.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/ring2nest.m (revision 25367) @@ -0,0 +1,51 @@ +function nPix = ring2nest(nSide,rPix) +% RING2NEST Convert MEALPix pixel numbers from ring to nest indexing +% +% nPix = ring2nest(nSide,rPix) +% +% +% nSide HEALPix resolution parameter +% rPix ring indexed MEALPix pixel numbers +% +% nPix next indexed MEALPix pixel numbers +% +% See also nest2ring +% +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: ring2nest.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + + +%% Parse and validate input +error(nargchk(1,2,nargin,'struct')); +if ~qNSide(nSide) + msgid = [mfilename ':nSide']; + error(msgid,msgid); +end +if isempty(rPix) + rPix = 1:(12*nSide^2); +elseif any(fix(abs(rPix(:))) ~= rPix(:)) || any(rPix(:)<1) + msgid = [mfilename ':rPix']; + error(msgid,msgid); +end + +% Convert +sz = size(rPix); +rPix = rPix - 1; +nPix = f_ring2nest(nSide,reshape(rPix,[],1)); +nPix = reshape(nPix+1,sz); + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/ring2z.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/ring2z.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/ring2z.m (revision 25367) @@ -0,0 +1,69 @@ +function z = ring2z(nSide,nRing) +% RING2Z Find HEALPix ring number corresponding to a given z coordinate +% +% z = ring2z(nSide,nRing) +% +% nSide HEALPix base-pixel sub-divisions +% nRing numeric array of ring numbers +% +% z size(nRing) numeric array of pixel center z coordinates +% +% Example +% z = ring2z(4,reshape(1:15,3,5)); +% +% See also z2ring +% +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: ring2z.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + + +%% Parse and validate input +p = inputParser; +p.addRequired('nSide',@qNSide); +p.addRequired('nRing',@(n)(qNRing(nSide,nRing))); +p.parse(nSide,nRing); + +%% pre-allocate result +z = zeros(size(nRing)); + +%% north polar cap +mask = nRing < nSide; +z(mask) = 1 - nRing(mask).^2/(3*nSide^2); + +%% tropical band +mask = ~mask & nRing < 3*nSide; +z(mask) = (2*nSide-nRing(mask))*2/(3*nSide); + +%% south polar cap +mask = nRing >= 3*nSide; +z(mask) = -1 + (4*nSide-nRing(mask)).^2/(3*nSide^2); + +return + +function tf = qNRing(nSide,nRing) +% QNRING true if valid ring number + +tf = isnumeric(nRing); +if tf + nRing = nRing(:); + tf = ... + all(nRing == fix(abs(nRing))) && ... + all(nRing > 0) && ... + all(nRing < 4*nSide); +end + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/vec2ang.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/vec2ang.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/vec2ang.m (revision 25367) @@ -0,0 +1,54 @@ +function tp = vec2ang(xyz) +% VEC2ANG Convert cartesian direction vector to position angle on sphere +% +% tp = vec2ang(xyz) +% +% xyz cell array of [3,1] cartesian coordinate direction vectors +% +% The direction vectors need not be normalized. +% +% tp cell array of [2,1] position angles ([theta; phi], radians) on +% sphere +% +% See also ang2vec +% +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: vec2ang.m 5831 2011-02-25 23:55:19Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + + +if ~iscell(xyz) + msgid = [mfilename ':xyz']; + error(msgid,msgid); +end +tp = cellfun(@v2a,xyz,'UniformOutput',false); + +return + +function tp = v2a(xyz) +% dVec2Ang un-nomralized cartesian vector to position angle on sphere +% +% tp = v2a(xyz) +% +% xyz [3,1] +% +% tp [2,1] of [theta, phi] + +assert(numel(xyz)==3,[mfilename ':v2a'],[mfilename 'v2a']); +xyz = xyz/norm(xyz); +tp = [acos(xyz(3)); atan2(xyz(2),xyz(1))]; + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/vec2pix.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/vec2pix.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/vec2pix.m (revision 25367) @@ -0,0 +1,147 @@ +function nPix = vec2pix(nSide,xyz,varargin) +% VEC2PIX Convert cartesian direction vectors to MEALPix pixel numbers +% +% nPix = vec2pix(nSide,xyz,'Param1',Value1,...); +% +% nSide HEALPix resolution parameter +% xyz cell array of [3,1] cartesian direction vectors +% +% Param Value +% 'nest' use nested indexing (true | {false}) +% +% nPix size(xyz) pixel number array +% +% Example +% xyz = {[1;0;0],[0;1;0],[0;0;1]}; +% nPix = vec2pix(4,xyz); +% +% See also ang2pix, pix2ang, pix2vec +% +% Requires ring2nest +% +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: vec2pix.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + + +%% Validate and parse input +p = inputParser; + +p.addRequired('nSide', ... + @(x)(isscalar(x) && isnumeric(x) && x>0 && nextpow2(x) == log2(x))); +p.addRequired('xyz',@iscell); +p.addParamValue('nest',false,@(n)(isscalar(n) && islogical(n))); + +p.parse(nSide,xyz,varargin{:}); +nest = p.Results.nest; + +%% Convert vector(s) to pixel numbers +% Process xyz as a vector then restore shape + +% Save xyz shape +sXyz = size(xyz); + +% Insure cell vectors are [3,1] +xyz = cellfun(@(v)(reshape(v,[],1)),xyz,'UniformOutput',false); + +% Re-express xyz as a column vector +xyz = cell2mat(reshape(xyz,1,[])); + +% Find pixels +nPix = f_vec2pix_ring(nSide,xyz(1,:),xyz(2,:),xyz(3,:)); +nPix = nPix + 1; +if nest + nPix = ring2nest(nSide,nPix); +end + +% Restore shape +nPix = reshape(nPix,sXyz); + +return + +function iPix = f_vec2pix_ring(nSide, x, y, zin) + +%======================================================================= +% renders the pixel number iPix (RING scheme) for a pixel containing +% a point on a sphere at coordinate vector (=x,y,z), given the map +% resolution parameter nside +%======================================================================= + +dnorm = arrayfun(@(x,y,z)(norm([x,y,z])),x,y,zin); +z = zin./dnorm; + +phi = atan2(y,x); +za = abs(z); +% phi(phi < 0.0) = phi(phi < 0.0) + 2*pi; +phi = mod(phi,2*pi); +tt = phi/(pi/2); +nl4 = fix(4*nSide); + +%% preallocate result +iPix = zeros(size(x)); + +%% Equatorial region +% Identify pixels +pixMask = 3*za <= 2; +if ~isempty(pixMask) + ttPix = tt(pixMask); + zPix = z(pixMask); + + temp1 = nSide*(ttPix + 0.5); + temp2 = nSide*zPix*0.75; + + % find ascending, descending line indices + jp = fix(temp1-temp2); + jm = fix(temp1+temp2); + + ir = fix(nSide+jp-jm); + kShift = fix(bitand(ir,1)); + ip = fix((jp+jm-nSide+kShift+1)/2); + ip(ip>=nl4) = fix(ip(ip>=nl4)-nl4); % can this be replaced by a mod? + iPix(pixMask) = fix(2*nSide*(nSide-1) + nl4*fix(ir) + ip); +end +%% North, south polar cap +% Identify pixels +pixMask = ~pixMask; +if ~isempty(pixMask) + ttPix = tt(pixMask); + tp = ttPix - fix(ttPix); + tmp = nSide*sqrt(3*(1-za(pixMask))); + + % find increasing, decreasing edge line indices + jp = fix(tp.*tmp); + jm = fix((1-tp).*tmp); + + % ring number counted from the closest pole + ir = fix(jp+jm+1); + ip = fix(ttPix.*ir); + tMask = ip>=4*ir; + ip(tMask) = fix(ip(tMask)-4*ir(tMask)); + + % north pole pixels + nPixMask = pixMask & z > 0; + tMask = z(pixMask) > 0; + iPix(nPixMask) = fix(2*ir(tMask).*(ir(tMask)-1) + ip(tMask)); + + % south pole pixels + sPixMask = pixMask & z <= 0; + tMask = ~tMask; + iPix(sPixMask) = fix(... + 3*nSide*fix(nl4) - 2*ir(tMask).*(ir(tMask)+1)+ip(tMask)... + ); +end + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/xy2pix.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/xy2pix.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/xy2pix.m (revision 25367) @@ -0,0 +1,72 @@ +function ipix = xy2pix(nSide,ix,iy,nFace) +% XY2PIX gives the pixel number ipix (NESTED) of ix, iy and face_num +% +% ipix = xy2pix(nSide, ix, iy, nFace) +% +% nside HEALPix resolution parameter (power of 2) +% ix,iy x,y location on the Healpix face (1 <= ix, iy <= nSide) +% nFace HEALPix face number (1 <= nFace <= 12) +% +% ipix MEALPix NESTED pixel index for given x,y on the face +% +% See also pix2xy +% +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: xy2pix.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +%% Parse and validate arguments +p = inputParser; +p.KeepUnmatched = false; + +p.addRequired('nSide',@qNSide); +qNXY = @(n)(all(n(:)>0) && all(n(:)<=nSide)); +p.addRequired('ix',qNXY); +p.addRequired('iy',qNXY); +p.addRequired('nFace',... + @(n)(isscalar(n) && isnumeric(n) && ... + n == fix(n) && 0 < n && n <= 12)); +p.parse(nSide,ix,iy,nFace); + +%% Initialization +[x2pix y2pix]=mk_xy2pix; + +%% Convert to HEALPix indexing +ix = ix - 1; +iy = iy - 1; + +scale = 1; +scaleFactor = 16384; +ipf = 0; +ismax = 1; +if nSide > 1048576 + ismax = 3; +end +for k = 1:ismax + ix_low = mod(ix,128); + iy_low = mod(iy,128); + ipf = ipf + (x2pix(ix_low+1) + y2pix(iy_low+1))*scale; + scale = scale*scaleFactor; + ix = fix(ix/128); + iy = fix(iy/128); +end +ipf = ipf + (x2pix(ix+1) + y2pix(iy+1))*scale; + +ipix = ipf + nFace*nSide^2; + +ipix = ipix+1; % MEALPix numbering + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/ylm.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/ylm.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/ylm.m (revision 25367) @@ -0,0 +1,109 @@ +function v = ylm(nSide,L,M,varargin) +% YLM returns spherical harmonic basis Y_L^M on the HEALPix sphere +% +% v = ylm(nSide,L,M,'Param1',Value1,'Param2',Value2,...); +% +% nSide HEALPix resolution parameter (power of 2) +% L spherical harmonic degree (0 <= L) +% M spherical harmonic order (-M <= L <= M) +% +% Param Value +% 'real' return real values (true | {false}) +% 'nest' nested indexing flag (true | {false}) +% +% v spherical harmonics evaluated at pixel centers +% +% Example: +% +% % Plot Ylm for (L,M) = (4,2) on nSide = 16 HEALPix sphere +% pix = ylm(16,4,2); +% hp3d(pix); +% +% Requires pix2ang +% Required by alm2pix, pix2alm +% See also hp3d +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: ylm.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +%% Parse and validate input +p = inputParser; +p.KeepUnmatched = false; +p.FunctionName = mfilename; + +p.addRequired('nSide',@qNSide); +p.addRequired('L',@(x)(isnumeric(x) && round(x) == x && 0 <= x)); +p.addRequired('M',@(x)(isnumeric(x) && round(x) == x)); +p.addParamValue('real',false,@(x)(isscalar(x) && islogical(x))); +p.addParamValue('nest',false,@(x)(isscalar(x) && islogical(x))); + +p.parse(nSide,L,M,varargin{:}); +nSide = p.Results.nSide; +L = p.Results.L; +M = p.Results.M; +real = p.Results.real; +nest = p.Results.nest; + +% Check m is in range +assert(abs(M)<=L,[mfilename ':M'],[mfilename ':M']); + +% Setup HEALPix pixels to solve on +tp = pix2ang(nSide,'nest',nest); + +% Evaluate for positive M, using symmetry to find value for negative M +signM = sign(M); +M = abs(M); + +% Exploit HEALPix rings for legendre function calc +th = cellfun(@(tp)(tp(1)),tp); +[b, ~, n] = unique(th); +ringLmn = legendre(L,cos(b)); +if L~=0 + ringLmn = squeeze(ringLmn(M+1,:,:)); +end +Lmn = ringLmn(n); + +% Calculate normalization constant +%a1 = (2*L+1); % 4pi normalization +a1 = (2*L+1)/(4*pi); % orthonormalization (default) + +% a2 = factorial(L-M)/factorial(L+M); +a2 = 1/prod((L-M+1):(L+M)); + +C = sqrt(a1*a2); + +% Calculate complex Ymn +phi = cellfun(@(x)(x(2)),tp); +Ymn = C*Lmn.*exp(1i*M*phi); + +% Use symmetry to find value for negative M +if 0 > signM + Ymn = (-1)^M*conj(Ymn); +end + +% Set v +v = Ymn; + +% Convert to real form if true == real +if real && 0 ~= M + if 0 < signM + v = (Ymn+conj(Ymn))/sqrt(2); + else % 0 > signM + v = (-1)^M*(conj(Ymn)-Ymn)/sqrt(-2); + end +end + +return Index: /issm/trunk-jpl/externalpackages/mealpix/src/z2ring.m =================================================================== --- /issm/trunk-jpl/externalpackages/mealpix/src/z2ring.m (revision 25367) +++ /issm/trunk-jpl/externalpackages/mealpix/src/z2ring.m (revision 25367) @@ -0,0 +1,55 @@ +function nRing = z2ring(nSide,z) +%% Z2RING Find HEALPix ring numbers from direction vector z coordinate +% +% nRing = z2ring(nSide,z) +% +% nSide HEALPix resolution parameter +% z numeric array of z values +% +% nRing size(z) numeric array of ring numbers +% +% Description +% HEALPix pixels centers are arranged in 4*nSide-1 rings of constant z. +% Each ring spans a given range of z. z2ring returns the ring numbers that +% cover the input z values. +% +% Example: +% nRing = reshape(1:15,3,5); +% z = ring2z(4,nRing); +% +% See also ring2z +% +% Authors: Lee Samuel Finn, Matthew Kinsey +% Copyright 2010-2011 Lee Samuel Finn + +% $Id: z2ring.m 5795 2011-02-20 04:01:32Z lsf@GRAVITY.PSU.EDU $ + +% Copyright 2010-2011 by Lee Samuel Finn. +% This program file is part of MEALPix. It is licensed under the Apache +% License, Version 2.0 (the "License"); you may not use MEALPix except in +% compliance with the License. You may obtain a copy of the License at +% +% +% Unless required by applicable law or agreed to in writing, MEALPix +% software distributed under the License is distributed on an "AS IS" +% BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +% implied. See the License for the specific language governing permissions +% and limitations under the License. + +%% Parse and validate input +p = inputParser; +p.FunctionName = mfilename; +p.addRequired('nSide',@qNSide); +p.addRequired('z',@(x)(all(abs(x(:))<=1))); +p.parse(nSide,z); + +%% Work +% Save input size +% Flatten input +% Find ring numbers +% Reshape ring numbers to reflect input + +zSz = size(z); +nRing = reshape(ringNum(nSide, z(:)),zSz); + +return