Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Just some minor performance improvements and warning fixes #2

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Just some minor performance improvements and warning fixes #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 6 additions & 4 deletions addmon.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,9 +46,11 @@
% [dems,dels,mz,lmcosi,mzi,mzo,bigm,bigl,rinm,ronm,demin]=addmon(L,m);
% dems(demin)
%
% Last modified by fjsimons-at-alum.mit.edu, 1/20/2011

defval('m',NaN)
% Previously modified by fjsimons-at-alum.mit.edu, 1/20/2011
% Minor changes by pfaff-at-kit.edu, 10/04/2016
if nargin==1
m=NaN;
end

matr=(repmat(0:L,2,1)'*diag([0 1]))';
dems=matranges(matr(:)')';
Expand Down Expand Up @@ -82,7 +84,7 @@
if L>=2
% Must do the first couple of degrees myself since gamini(x,0) is bad
lopos=cumsum(2*els(1:end-1)+1)+1;
ka=[repmat(1,1,L-1) ; els(3:end)-1 ; repmat(1,1,L-1); els(3:end)];
ka=[ones(1,L-1) ; els(3:end)-1 ; ones(1,L-1); els(3:end)];
ko=[[1 0 -1 2] gamini(repmat([0 -2 -1 2],1,L-1),ka(:)')]';
ko(lopos)=[2:2:2*L];
rinm=cumsum(ko);
Expand Down
12 changes: 8 additions & 4 deletions addmup.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,8 @@
%
% See also ADDMOFF, ADDMOUT
%
% Last modified by fjsimons-at-alum.mit.edu, 05/17/2011
% Previously modified by fjsimons-at-alum.mit.edu, 05/17/2011
% Last modified by Florian Pfaff for libDirectional, 10/04/2016

% Calculates the number of real spherical harmonic orders
% that belong to an expansion from degree l=0 to L - or vice versa.
Expand All @@ -34,9 +35,12 @@
% scalars. Note: degrees 0 and 1 are of course included.

% Check the behavior for the negatives should you ever need those

defval('nr',3)
defval('drk','a')
if nargin==0
nr=3;
drk='a';
elseif nargin==1
drk='a';
end

switch drk
case 'a'
Expand Down
12 changes: 7 additions & 5 deletions gamini.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,12 +15,14 @@
% See also DEGAMINI,GAMINI2
%
% Last modified by fjsimons-at-alum.mit.edu, 03/28/2009

defval('folding',3)
% Minor changes by pfaff-at-kit.edu, 10/04/2016
if nargin==1
folding=3;
end

% Copy a single scalar folding for all elements in the array
if prod(size(folding))==1
folding=repmat(folding,prod(size(data)),1);
if numel(folding)==1
folding=repmat(folding,numel(data),1);
end

% Never had a replication factor of 0 before
Expand All @@ -35,7 +37,7 @@
folding=folding(:)';

if ~all(size(data)==size(folding))
error([ ' Sizes of input and folding must be the same'])
error('Sizes of input and folding must be the same')
end

gelp=zeros(1,sum(folding));
Expand Down
81 changes: 55 additions & 26 deletions xyz2plm.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,19 +54,38 @@
%
% See also PLM2XYZ, PLM2SPEC, PLOTPLM, etc.
%
% Last modified by fjsimons-at-alum.mit.edu, 09/04/2014
% Previously modified by fjsimons-at-alum.mit.edu, 09/04/2014
% Minor changes by pfaff-at-kit.edu, 10/04/2016

t0=clock;

defval('method','im')
defval('lon',[])
defval('lat',[])
defval('dw',[])
defval('cnd',[])
switch nargin
case 1
method='im';
lat=[];
lon=[];
cnd=[];
case 2
method='im';
lat=[];
lon=[];
cnd=[];
case 3
lat=[];
lon=[];
cnd=[];
case 4
lon=[];
cnd=[];
case 5
cnd=[];
end

dw=[];

as=0;
% If no grid is specified, assumes equal spacing and complete grid
if isempty(lat) & isempty(lon)
if isempty(lat) && isempty(lon)
% Test if data is 2D, and periodic over longitude
fthph=reduntest(fthph);
polestest(fthph)
Expand Down Expand Up @@ -116,12 +135,14 @@
end

% Decide on the Nyquist frequency
defval('L',Lnyq);
if nargin==1
L=Lnyq;
end
% Never use Libbrecht algorithm... found out it wasn't that good
defval('libb',0)
libb=0;
%disp(sprintf('Lnyq= %i ; expansion out to degree L= %i',Lnyq,L))

if L>Lnyq | nlat<(L+1)
if L>Lnyq || nlat<(L+1)
warning('XYZ2PLM: Function undersampled. Aliasing will occur.')
end

Expand All @@ -143,14 +164,23 @@
error('Specify valid method')
end

fnpl=sprintf('%s/LSSM-%i-%i.mat',...
fullfile(getenv('IFILES'),'LEGENDRE'),L,length(x));

if exist(fnpl,'file')==2 & as==1

% Set to true to use the folder in which xyz2plm resides as the folder for the LEGENDRE folder
createInFolderOfMatlabFile=false;
if createInFolderOfMatlabFile
mfn=mfilename('fullpath'); %#ok<UNRCH>
fnpl=fullfile(mfn(1:end-8),'LEGENDRE',sprintf('LSSM-%i-%i.mat',L,length(x)));
else
fnpl=sprintf('%s/LSSM-%i-%i.mat',...
fullfile(getenv('IFILES'),'LEGENDRE'),L,length(x));
end


if exist(fnpl,'file')==2 && as==1
load(fnpl)
else
% Evaluate Legendre polynomials at selected points
Plm=repmat(NaN,length(x),addmup(L));
Plm=NaN(length(x),addmup(L));
if L>200
h=waitbar(0,'Evaluating all Legendre polynomials');
end
Expand Down Expand Up @@ -180,14 +210,14 @@
case {'irr'}
Plm=[Plm.*cos(lon(:)*m(:)') Plm.*sin(lon(:)*m(:)')];
% Add these into the sensitivity matrix
[C,merr,mcov,chi2,L2err,rnk,dw]=datafit(Plm,fthph,[],[],cnd);
[C,merr,mcov,chi2,L2err,rnk,dw]=datafit(Plm,fthph);
lmcosi(:,3)=C(1:end/2);
lmcosi(:,4)=C(end/2+1:end);
case {'im','gl','simpson'}
% Perhaps demean the data for Fourier transform
defval('dem',0)
dem=false;
if dem
meanm=mean(fthph,2);
meanm=mean(fthph,2); %#ok<UNRCH>
fthph=fthph-repmat(meanm,1,nlon);
end

Expand All @@ -200,7 +230,7 @@

if dem
% Add the azimuthal mean back in there
gfft(:,1)=2*pi*meanm;
gfft(:,1)=2*pi*meanm; %#ok<UNRCH>
end

% Note these things are only half unique - the maximum m is nlon/2
Expand All @@ -220,8 +250,8 @@
a(:,1)=a(:,1)/2;
b(:,1)=b(:,1)/2;
Pm=Plm(:,mz(ord+1:end)+ord)*pi;
[lmcosi(mz(ord+1:end)+ord,3)]=datafit(Pm,a(:,ord+1),[],[],cnd);
[lmcosi(mz(ord+1:end)+ord,4)]=datafit(Pm,b(:,ord+1),[],[],cnd);
[lmcosi(mz(ord+1:end)+ord,3)]=datafit(Pm,a(:,ord+1));
[lmcosi(mz(ord+1:end)+ord,4)]=datafit(Pm,b(:,ord+1));
end
case 'simpson'
% Loop over the degrees. Could go up to l=nlon if you want
Expand Down Expand Up @@ -250,7 +280,6 @@
lmcosi(addmup(l-1)+1:addmup(l),3)=clm(:)/4/pi;
lmcosi(addmup(l-1)+1:addmup(l),4)=slm(:)/4/pi;
end
rnk=[];
end

% Get rid of machine precision error
Expand All @@ -264,8 +293,8 @@
% Tests if last longitude repeats last (0,360)
% and removes last data column
if sum(abs(grd(:,1)-grd(:,end))) >= size(grd,2)*eps*10
disp(sprintf('Data violate wrap-around by %8.4e',...
sum(abs(grd(:,1)-grd(:,end)))))
fprintf('Data violate wrap-around by %8.4e\n',...
sum(abs(grd(:,1)-grd(:,end))));
end
grd=grd(:,1:end-1);

Expand All @@ -274,6 +303,6 @@ function polestest(grd)
% Tests if poles (-90,90) are identical over longitudes
var1=var(grd(1,:));
var2=var(grd(end,:));
if var1>eps*10 | var2>eps*10
disp(sprintf('Poles violated by %8.4e and %8.4e',var1,var2))
if var1>eps*10 || var2>eps*10
fprintf('Poles violated by %8.4e and %8.4e\n',var1,var2);
end