{\rtf1\ansi\ansicpg1252\deff0\deflang1033\deflangfe1033{\fonttbl{\f0\fswiss\fprq2\fcharset0 Helvetica;}{\f1\froman\fprq2\fcharset0 Times New Roman;}}
{\*\generator Msftedit 5.41.15.1503;}\viewkind4\uc1\pard\nowidctlpar\tx720\tx1440\tx2160\tx2880\tx3600\tx4320\tx5040\tx5760\tx6480\tx7200\tx7920\tx8640\f0\fs24 function shdemo(l,m)\par
% Type\par
% SHDEMO(l,m) and see the cos-Plm plotted on a sphere\par
%\par
% \par
% \par
\par
if  (m > l | m <0)\par
  error('Incorrect degree/order specification')\par
end\par
\par
cosorsin=1;  %2 for sine\par
bg=0;\par
ampl=1;\par
\par
coeffs=repmat(bg,addmup(l,'a'),2);\par
coeffs(addmup(l-1,'a')+1+m,cosorsin)=ampl;\par
coeffs=coeffs';\par
\par
% Adaptation for Oded's program\par
data=sha2grid(coeffs(:));\par
\par
% Do something with the output\par
plotplm(data);\par
\par
\par
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\par
function plotplm(data)\par
% PLOTPLM(data)\par
\par
% \par
\par
% Make sphere for the data\par
lons=linspace(0,360,100)/180*pi;\par
lats=linspace(90,-90,100)/180*pi;\par
[lons,lats]=meshgrid(lons,lats);\par
rads=ones(size(lats));\par
[x,y,z]=sph2cart(lons,lats,rads);\par
\par
surface(x,y,z,'FaceColor','texture','Cdata',data);\par
\par
axis image\par
shading flat\par
view(220,25)\par
\par
set(gca,'color',[.9 .9 .9]-.1)\par
axis('off')\par
colormap(gray(2))\par
\par
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\par
function ditist=addmup(n,allertour)\par
% ditist=ADDMUP(number,direction)\par
%\par
% Where direction is 'a' or 'r'\par
% 'a' calculates the number from input l\par
% 'r' calculates the l from input number\par
%\par
% m=0:l therefore the number m is l+1\par
% This program calculates \\sum\\limits_\{l=0\}^\{L\}(l+1) for a given L.\par
%\par
% This function calcuates the number of real spherical harmonic order m\par
% that belong to a certain degree number l - or vice versa, depending on the \par
% specification 'a' or 'r'. Useful and used in PLMSPEC.\par
% Works for matrices, vectors as well as scalars.\par
% Note how we very explicitly include degree 0 and 1 - this may\par
% have to be remediated in function calls that work in a coordinate \par
% frame of which the origin coincides with the center of mass!\par
\par
% Last modified September 27th 2000\par
\par
switch  allertour\par
  case 'a'\par
    ditist=1/2*(n+1).^2+1/2*n+1/2;\par
  case 'r'\par
    ditist=-3/2+1/2*sqrt(1+8*n);\par
    if prod(round(ditist)==ditist)==0\par
      error('Invalid entry - noninteger result')\par
    end\par
end\par
\par
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\par
function varargout=sha2grid(clmi, lmin, lmax, res, unit)\par
% sha2grid: function to convert spherical harmonic coefficients to spatial grid.\par
%\par
% Calling sequence:\par
%  [r,lon,lat]=sha2grid(clmi, lmin, lmax, res[, unit]);\par
%  or \par
% r=sha2grid(clmi, lmin, lmax, res[, unit]);\par
% \par
% clmi is an array of coefficients clmi=[A_l^m B_l^m]; clmi=clmi'(:) for all l and m\par
% lmax is the highest degree desired in the expansion, which defaults to the\par
% highest degree supplied. lmin is the minimum l desired which defaults to 0.\par
% res is the resolution of the output grid in degrees, which defaults to 1.\par
% unit is the units of the output lat,lon arrays ('deg' or 'rad'),  which \par
% defaults to 'deg'.\par
\par
% Oded Aharonson; 5/1/00\par
\par
defaultval('unit','deg');\par
defaultval('res',1);\par
defaultval('lmax',ind2lmi(length(clmi)));\par
defaultval('lmin',0);\par
\par
lon=linspace(0,360,360/res+1);\par
lat=linspace(-90,90,180/res+1);\par
\par
lon=lon*pi/180; lat=lat*pi/180;\par
nlon=length(lon); nlat=length(lat);\par
\par
\lang1031 grid=zeros(nlat,nlon);\par
\par
x=sin(lat);\par
\par
\lang1033 for l=lmin:lmax,\par
  ylm=lgndr(l,x);\par
  m=0:l;\par
  clm=clmi(lmi2ind(l,m,1));\par
  slm=clmi(lmi2ind(l,m,2));\par
\par
  mlon=mtimes(m',lon);\par
\par
  fac1=repmat(clm,1,nlon).*cos(mlon)+repmat(slm,1,nlon).*sin(mlon);\par
  \par
  grid=grid+ylm'*fac1;\par
end\par
\par
varargout\{1\}=grid;\par
if nargout == 3,\par
  if strcmp(unit,'deg')\par
    lon=lon*180/pi;\par
    lat=lat*180/pi;\par
  elseif strcmp(unit,'rad')\par
    % do nothing\par
  else\par
    disp('unit argument must be "deg" or "rad". Assumed "rad"');\par
  end\par
  [varargout\{2\},varargout\{3\}]=meshgrid(lon,lat);\par
end\par
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\par
function defaultval(arg,defval,quiet)\par
% defaultval: sets a variable to a default value if the variable is undefined.\par
%\par
% Checks if a variable arg exists in the caller's workspace,\par
% if nonexistent or if empty arg is set to value defval in caller's workspace\par
% If quiet does not exist or is 1, then will print value being used to screen.\par
\par
% Written by Keli, 5/1/00\par
% Modified by Oded Aharonson, 5/1/00\par
\par
if ~exist('quiet'),\par
  quiet=0;\par
end\par
\par
flag=1;\par
if evalin('caller',[ 'exist(''' arg ''')']);,\par
  flag=evalin('caller',[ 'isempty(' arg ')']);\par
end\par
if flag,\par
  if ~quiet,\par
    if (ischar(defval) | length(defval)<3)\par
      [path,prog,ext,ver]=star69;\par
      verbose([prog,': Using ',arg,'=',num2str(defval)]);\par
    end\par
  end\par
  assignin('caller',arg,defval);\par
end\par
  \par
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\par
function [path,name,ext,ver]=star69()\par
% star69: who called me? returns calling function name if exists.\par
\par
str=dbstack;\par
if length(str) > 2\par
  str=str(3).name;\par
else\par
  str='';\par
end\par
[path,name,ext,ver]=fileparts(str);\par
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\par
function varargout=verbose(str,level)\par
% verbose: function to display message str if verbose is turned on.\par
%\par
% By default verbose is on. To turn off use\par
% verbose off;\par
% verbose with no arguments returns the current verbose state.\par
% Advance usage allows verbose levels: off is 0, on is 1, debug is 2,\par
% and higer debug modes follow.\par
%\par
% Examples:\par
% verbose on;         % turns verbose on  (level=1)\par
% verbose off;        % turns verbose off (level=0)\par
% verbose(str);       % simple usage.\par
% verbose(str,0);     % print a "level 0" message.\par
% verbose('debug')    % set verbose level to 2.\par
% verbose(['x: ',num2str(x)],'debug')   % Print a debug message.\par
% verbose             % return current verbose level.\par
\lang1031 %\par
% Oded Aharonson; 5/3/00\par
\par
persistent verbose_flag ;\par
\lang1033 %mlock;\par
\par
if (length(verbose_flag) == 0)\par
  verbose_flag=1;\par
end\par
\par
if (exist('level') & ischar(level))\par
  if strcmp(level,'off')\par
    level=0;\par
  elseif strcmp(level,'on')\par
    level=1;\par
  elseif strcmp(level,'debug')\par
    level=2;\par
  else\par
    error('verbose: illegal label!')\par
  end\par
end\par
  \par
if ~exist('str'),\par
  varargout\{1\}=verbose_flag;\par
\lang1031 else\par
  if ~ischar(str),\par
    verbose_flag=str;\par
  \lang1033 elseif strcmp(str,'off'),\par
    if ~exist('level'),\par
      \lang1031 verbose_flag=0;\par
    else\par
      if level<=verbose_flag,\par
\tab\lang1033 disp(str);\par
      end\par
    end\par
  elseif strcmp(str,'on'),\par
    if ~exist('level'),\par
      verbose_flag=1;\par
    else\par
      if level<=verbose_flag,\par
\tab disp(str);\par
      end\par
    end\par
  elseif strcmp(str,'debug'),\par
    if ~exist('level'),\par
      verbose_flag=2;\par
    else\par
      if level<=verbose_flag,\par
\tab disp(str);\par
      end\par
    end\par
  else\par
    if ~exist('level'), level=1; end\par
    if level<=verbose_flag,\par
      disp(str);\par
    end\par
  end\par
end\par
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\par
\par
function [l,m,i]=ind2lmi(ind)\par
% ind2lmi: function to convert spherical harmonic (l,m,i) to index\par
% i=1,2, even for m=0\par
\par
if any(ind < 1),\par
  error('ind2lmi: index must be > 0!');\par
end\par
i=mod(ind-1,2)+1;\par
tmp=ind-i;\par
l=fix((-1+sqrt(1+4.*tmp))./2);\par
m=(ind-i-l.*(l+1))/2;\par
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\par
function y=lgndr(l,x)\par
% lgndr: function to return normalized legendre polymomials.\par
%\par
% Normalized such that:\par
% int_theta,phi lgndr(l)*cos(m*phi) sin(theta) dtheta dphi=4 pi\par
%\par
% Normalization Test:\par
%\par
% dth=.001; dphi=.001;\par
% l=2; m=1;\par
% th=0:dth:pi; phi=0:dphi:(m*pi);\par
% pl=lgndr(l,cos(th)); pl=pl(m+1,:);\par
% sum(pl.^2.*sin(th)*dth)*sum(cos(m*phi).^2*dphi)/(4*pi)\par
% \par
% Oded Aharonson, 05/15/00\par
% Modified 6/12/00, for better accuracy at large l.\par
\par
% For better accuracy at l>80 :\par
if size(x,2)==1,\par
  x=x';\par
end\par
m=0:l;\par
y=legendre(l,x,'sch');\par
y=y*sqrt(2*l+1);\par
\par
% Modified from:\par
% for i=1:length(m)\par
%  y(i,:) = y(i,:)*(-1)^(m(i))*sqrt( (2-(m(i)==0))*(2*l+1)*factorial(l-m(i)) / factorial(l+m(i)) );\par
%end\par
\par
\par
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\par
function ind=lmi2ind(l,m,i)\par
% lmi2ind: function to convert spherical harmonic (l,m,i) to index\par
% i=1,2, even for m=0\par
\par
if any((i<1) | (i>2)), error('lmi2ind: bad value for i!'); end\par
if any(m<0), error('lmi2ind: bad value for m!'); end\par
if any(l<0), error('lmi2ind: bad value for l!'); end\par
\par
la=[length(l),length(m),length(i)];\par
\par
if ( (prod(la) > max(la)) & any(abs(diff(la))>0) )\par
  [l,m,i]=ndgrid(l,m,i);\par
  l=l(:); m=m(:); i=i(:);\par
  f=find(m<=l);\par
  l=l(f); m=m(f); i=i(f);\par
else\par
  if any(m>l), error('lmi2ind: bad value for m!'); end\par
end\par
\par
ind=sort(l.*(l+1)+2.*m+i);\f1\par
}