Upload New File

get_poles_via_aaa.m

+function [poles,res] = get_poles_via_aaa(a,min_poles,norm_fkoeffs,delta)
+% poles = get_poles_via_aaa(a,min_poles,A);
+%
+% Der Vektor a enthaelt die relevanten Laurentkoeffizienten
+% Output sind alle Pole der rationalen Approximation (mit mindestens min_poles
+% Polen), die in der Naehe der Einheitskreislinie liegen. 
+% Maximal length(a)-1 Pole werden bestimmt.
+% A ist die L^2-Norm des Fernfelds (bzw. die l2-Norm aller Fourierkoeffs).
+
+M = length(a);
+
+rho = 1.0 * (4*norm_fkoeffs/delta)^(1/M);
+tol = (1/2) * 4/(M/2) * sqrt( (1-rho^(-2*M-2))/(1-rho^(-2)) )/sqrt(2*pi) * delta/rho;
+
+t = (0:99)/100*2*pi;
+zeta = rho*exp(1i*t);
+
+f = zeros(size(zeta));
+f_reverse = f;
+for n=1:M    
+   f = f + a(n)*zeta.^(-n);
+   f_reverse = f_reverse + conj(a(end-n+1))*zeta.^(-n);
+end
+
+
+pole_count = floor(M/2)*2;
+
+poles = [];
+res = [];
+
+fhelp = f;
+[poles_help,res_help] = special_aaa_infty(zeta,fhelp,0,tol,min_poles,pole_count);
+poles = poles_help;  
+res = res_help;
+
+fhelp = f_reverse;
+[poles_help,res_help] = special_aaa_infty(zeta,fhelp,0,tol,min_poles,pole_count);
+poles = [poles; poles_help];
+res = [res; res_help];
+
+fhelp = (f+f_reverse)/2;
+[poles_help,res_help] = special_aaa_infty(zeta,fhelp,0,tol,min_poles,pole_count);
+poles = [poles; poles_help];
+res = [res; res_help];
+
+fhelp = (f-f_reverse)/2;
+[poles_help,res_help] = special_aaa_infty(zeta,fhelp,0,tol,min_poles,pole_count);
+poles = [poles; poles_help];
+res = [res; res_help];
+
+
+%% eliminiere alle Pole, deren Radius nicht im Intervall [0.95,1.05] liegen
+abs_poles = abs(poles);
+inds = find((abs_poles<0.95) | (abs_poles>1.05));
+poles(inds) = [];
+res(inds) = [];
+
+
+
+